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ABSTRACT 


A  constitutive  microdamage  model  is  developed  capable  of  simulating  high  shock 
compression,  release,  dilatation  (tension),  and  microdamage  evolution  leading  possibly  to 
fracture  and  penetration  of  targets  after  hypervelocity  impact.  The  microdamage 
constitutive  model  is  applicable  to  polycrystalline  metals  and  is  appropriate  in  the  lower 
range  of  hypervelocity  impact  velocity,  i.e.  approximately  2-7  Km/s,  over  which  the 
projectile  and  target  materials  remain  in  the  solid  state.  The  model  implements  the  Mie- 
Gruneisen  equation  of  state  coupled  with  the  Hugoniot  relations  along  with  expressions 
of  non-linear  elastic  moduli  (bulk  and  shear)  as  functions  of  volume  strain,  temperature 
and  microdamage.  The  viscoplastic  material  response  includes  strain  and  strain  rate 
hardening  and  temperature  and  microdamage  softening.  The  microdamage  evolution 
model  is  based  on  the  micromechanics  of  an  expanding  void,  and  is  capable  of  modeling 
void  compaction  and  expansion  that  leads  to  spall-fracture  as  an  evolutionary  time 
dependent  process.  The  constitutive  microdamage  model  was  implemented  in  the 
Autodyn™  software  and  a  series  of  computer  simulations  of  hypervelocity  impact 
experiments  on  Alnoo  plates  with  soda-lime  glass  spherical  projectiles  were  conducted. 
The  results  of  the  simulations  are  compared  with  the  laboratory  experimental  results  in 
terms  of  crater,  penetration  hole  and  back-wall  spallation  geometry  of  the  target  plate. 
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STATEMENT  OF  PROBLEM 


The  importance  of  a  reliable  constitutive  model  to  describe  the  deformation 
process  and  the  eventual  material  degradation  and  fracture  after  hypervelocity  impact  has 
become  of  primary  interest  in  the  last  decade.  The  deployment  of  low  earth  orbit 
satellites  capable  of  surviving  high  velocity  impacts  by  micro-meteoroids  and  small 
orbital  space  debris  has  been  and  continues  to  be  the  subject  of  extensive  research.  In 
addition,  the  eventual  construction  of  the  International  Space  Station  will  put  to  test  the 
extensive  experimental  and  theoretical  research  being  developed  in  the  area  of 
hypervelocity  impact  phenomena.  Furthermore,  the  importance  of  this  area  of  research  is 
not  limited  to  the  aerospace  industry.  The  US  Defense  Department  has  always 
recognized  the  relevance  of  an  in-depth  understanding  of  the  armor  penetration  process  to 
aid  in  the  design  of  armor  structures  and  armor  penetration  devices. 

In  recent  years  there  have  been  significant  research  efforts  directed  towards  the 

development  of  constitutive  models  that  simulate  crater  formation,  spall  fracture  and 

penetration  as  a  result  of  a  hypervelocity  impact.  Many  of  these  models  have  been 

implemented  and  built  into  several  public  domain  and  commercial  software  packages' that 

are  available  today  (Steinberg  et  al.  1980  and  1989,  Johnson  et  al.  1983,  Zerilli  et  al. 

1987).  However,  these  constitutive  models  have  known  limitations.  The  most  important 

limitation  has  to  do  with  the  modeling  of  the  dynamic  failure  process  where,  in  most 

cases,  fracture  is  assumed  to  take  place  instantaneously  without  prior  accumulation  of 

material  micro-degradation.  The  development  of  more  sophisticated  constitutive  models 

vii 


that  can  describe  failure  is  crucial  for  the  problem  of  modeling  the  physical  processes 
associated  with  hypervelocity  impact  events,  especially  in  relation  to  the  design  of  armor, 
armor  penetrators,  debris-shield  design  and  other  applications  where  the  lethality  of  the 
impact  is  a  design-controlling  factor  (e.g.,  hypervelocity  impact  destruction  of  incoming 
ballistic  warhead). 

It  is  the  goal  of  this  dissertation  to  develop  a  constitutive-microdamage  model  that 
can  describe  the  thermo-mechanical  deformation  associated  with  high  shock 
compression,  release,  dilatation  and  shear,  and  the  progressive  microdamage  that  leads  to 
fracture  of  ductile  polycrystalline  materials. 
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CHAPTER  1 


INTRODUCTION 


The  physical  features  of  projectile-target  collisions  at  very  high  velocities  are 
varied  and  depend  upon  the  interplay  of  major  controlling  geometric  and  physical 
parameters  such  as  the  velocity  of  impact,  the  relative  dimensions  of  the  target  thickness 
and  the  projectile  diameter,  and  the  material  properties  of  the  target  and  the  projectile. 
There  are  many  damage  scenarios  that  can  be  identified,  ranging  from  target  crater 
formation  to  target  perforation,  where  for  the  latter  situation  at  very  high  projectile 
velocities  it  is  possible  to  have  target  and  projectile  fragmentation  and  debris  cloud 
formation,  while  for  still  higher  velocities  melting  and  vaporization  of  portions  of  the 
projectile  and  target  materials  are  also  possible  (Anderson  et  al.  1993,  Horz  et  al.  1995). 

Spherical  projectiles  traveling  at  velocities  varying  from  2  to  7  km/s  can,  upon 
impact,  produce  spherical  shock  waves  of  compression  and  shear  with  pressures  ranging 
from  10-150  GPa,  that  propagate  into  both  the  target  and  projectile  at  very  high  strain 
rates  of  the  order  of  106-108  s'1.  Transmission  electron  microscopy  has  shown  evidence 
of  varied  arrays  of  deformation-induced  defects  and,  in  some  cases,  dynamic 
recrystallization  immediately  below  the  impact  surface  of  several  metallic  targets  at  the 
higher  end  of  the  impact  velocity  range,  indicating  substantial  temperature  increases  in 
both  the  target  and  projectile  (Rivas  et  al.  1995,  Quinones  et  al.  1998).  The  transit  time 
of  the  wave  propagation  and  thermomechanical  material  response  across  the  target 
thickness  is  of  the  order  a  couple  of  microseconds  (ps),  implying  essentially  adiabatic 


l 
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deformation.  Because  of  thermomechanical  coupling  associated  with  the  substantial 
dilatational  and  shear  deformations,  the  highly  deformed  thermally  softened  region 
beneath  the  contact  surface  of  the  target  provides  the  displaced  material  that  forms  the 
surface  crater  ring  and  ejecta.  For  the  target  material  below  the  crater,  a  competing  effect 
develops  between  the  strain  hardening  and  thermal  softening  effects  as  evidenced  by  the 
measured  vickers  microhardness  profiles  (Quinones  et  al.  1998).  Further  across  the 
target  thickness,  the  attenuated  shock  wave  produces  a  reflected  tensile  wave  traveling 
back  from  the  rear  surface  of  the  target.  If  the  tensile  mean  stress  reaches  a  threshold 
value  the  growth  and/or  nucleation  of  voids  can  begin,  producing  mechanical  softening  of 
the  material  and  possibly  coalescence  of  voids  leading  to  fracture  (Curran  et  al.  1987). 

The  development  of  computer  codes  capable  of  predicting  the  deformation  and 
damage  processes  produced  by  hypervelocity  impact  have  evolved  from  being  very 
rudimentary  in  the  50’s  and  60’s,  to  very  sophisticated  at  this  time.  One  of  the  initial 
limitations  of  these  codes  was  the  absence  of  shear  deformation  considerations  in  the 
simulation  algorithms.  The  development  of  more  sophisticated  instrumentation  for  the 
experimental  characterization  of  the  deformation  process  under  shock  loading  (Meyers 
1994)  has  improved  our  understanding  of  the  hypervelocity  impact  problem,  and  with  the 
advent  of  high  speed  computers  that  use  more  sophisticated  material  models  and 
numerical  algorithms,  the  applicability  of  the  computer  codes  has  been  extended.  Now, 
computer  simulations  of  high-rate  deformation  processes  are  being  used  in  industry  as 
well  as  in  military  and  aerospace  applications  as  part  of  the  design  process. 
Unfortunately,  one  of  the  limitations  that  still  lingers  in  high  shock  compression 
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computer  codes  is  the  poor  constitutive  modeling  of  the  dynamic  fracture  process,  where 
it  is  often  assumed  that  fracture  occurs  instantaneously  without  prior  evolution  of  the 
material  degradation. 

1.1  Background 

The  following  is  a  brief  discussion  of  the  research  articles  that  form  the  base  for 
the  research  work  presented  in  this  dissertation.  The  literature  review  is  presented  in 
chronological  order. 

1.1.1  Constitutive  Models 

The  articles  discussed  in  this  section  deal  specifically  with  two  topics:  pressure- 
volume  relations  and  yield  models  (in  the  Von  Mises  sense)  of  polycrystalline  metals. 
The  pressure-volume  relations  discussed  are  for  pressures  up  to  10  TPa.  The  dynamic 
yield  models  varied  in  complexity  from  perfectly  plastic  rate-independent  to  viscoplastic 
with  strain  hardening  and  thermal  softening. 

Walsh  et  al.  (1957)  presented  an  extensive  compilation  of  material  parameters 
that  describe  the  Hugoniot  curves,  Mie-Gruneisen  pressure  volume  relations  and  the 
volume  strain  dependence  of  the  Gruneisen  coefficient  for  shock  pressures  of  up  to  400 
kbars  (40GPa). 

Kratochvil  et  al.  (1969,1970)  proposed  an  elastic-viscoplastic  constitutive  model 
that  described  strain  and  strain  rate  hardening  as  thermodynamic  functions  that  described 
the  movement  and  interaction  of  microstructural  defects. 
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Davison  et  al.  (1979)  presented  a  very  thorough  exposition  of  the  then  status  of 
investigations  concerning  the  response  of  solids  to  shock  compression.  They  described 
the  mechanical  and  structural  behavior  of  materials  under  shock  loading  beginning  with 
the  elastic  response,  the  pressure-volume  relations,  hydrodynamic  approximation  to  the 
response  of  solids,  plastic  response  and  finally  spall  fracture.  They  also  described  the 
electrical,  magnetic  and  optical  properties  of  solids  under  shock  loading. 

Steinberg  et  al.  (1980)  proposed  a  constitutive  model  applicable  at  high  strain 
rates.  The  constitutive  equation  models  viscoplastic  deformation  including  descriptions 
of  the  shear  modulus  and  yield  stress  that  incorporate  dependence  on  pressure, 
temperature  and  plastic  strain.  However,  they  restricted  the  model  to  high  strain  rate 
(>105s  ')  deformation  by  assuming  a  rate  independent  model.  Steinberg  et  al.  (1988) 
developed  an  addition  to  the  Steinberg-Guinan  model  to  extend  its  validity  to  plastic 
strain  rates  as  low  as  HTV1. 

Godwal  et  al.  (1983)  presented  a  review  of  pressure-volume  relation  of  condensed 
matter  at  up  to  10  TPa.  They  described  many  pressure-volume  relations  for  a  wide  range 
of  densities  and  temperatures,  obtained  by  empirical  and  quantum  mechanical 
calculations,  where  a  variety  of  physical  phenomena  such  as  lattice  and  electronic 
thermal  excitation,  phase  transition,  etc.  are  encountered. 

Johnson  et  al.  (1983)  developed  an  empirical  yield  model  for  metals  under  high 
strain,  high  strain  rate  and  high  temperature.  They  also  presented  a  database  of  material 
parameters  fitted  to  their  model  based  on  torsion,  Hopkinson  bar  and  quasi-static  tensile 
experiments.  This  model,  known  as  the  Johnson-Cook  model,  is  widely  used  in  high- 
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shock  computer  codes.  Rule  et  al.  (1998)  presented  a  modification  to  the  Johnson-Cook 
model  to  account  for  the  sudden  jump  in  magnitude  of  the  strain  rate  effect  that  many 
metals  exhibit  at  rates  greater  than  10V.  Follansbee  et  al.  (1988)  concluded  that  the 
increased  strain  rate  sensitivity  -was  due  to  a  dramatic  increase  in  dislocation 
accumulation  at  strain  rates  exceeding  ~10V,  and  proposed  a  model  that  used  the 
dislocation  density  as  the  only  state  variable  to  describe  the  plastic  behavior  of  f.c.c. 
(Face  center  cubic)  metals. 

Doraivelu  et  al.  (1984)  proposed  a  yield  function  model  for  porous  materials.  In 
their  model  the  yield  function  was  explicitly  coupled  to  the  first  and  second  invariants  of 
the  stress  tensor. 

Asay  et  al.  (1987)  presented  a  very  brief  discussion  of  the  basic  theoretical 
concepts  that  form  the  basis  for  the  development  of  material  models.  They  also  discussed 
some  experimental  techniques  used  in  the  characterization  of  the  shock  response  of 
materials  and  a  basic  discussion  of  the  pressure-volume  response  of  materials  in  the  solid, 
liquid  and  gaseous  states.  Finally  they  gave  a  very  basic  discussion  of  the  viscoplastic 
behavior  of  materials  and  some  fracture  considerations. 

Zerilli  et  al.  (1987)  presented  a  viscoplastic  constitutive  model  that  is  based  on 
dislocation  mechanics.  Their  model  applies  to  fee  poly  crystal  line  metals  and 
incorporates  the  effects  of  thermally  activated  mechanisms  of  dislocation  movement  and 
the  influence  of  solute  and  grain  size.  They  also  presented  material  parameters  for  OFHC 
Copper  and  Armco  Iron  based  on  the  data  presented  by  Johnson  et  al.  (1983). 

Wang  et  al.  (1997)  presented  a  study  of  the  applicability  of  three  viscoplastic 
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models  to  the  modeling  of  instabilities  due  to  strain  and  strain  rate  softening.  They 
considered  the  Perzyna,  Duvaut-Lions  models  and  a  new  model  proposed  by  them.  They 
also  presented  a  finite  difference  formulation  of  the  models  for  inclusion  in  an  implicit  or 
explicit  finite  element  code. 

Fossum  (1998)  presented  an  optimization  technique  to  fit  experimental  data  to  a 
viscoplastic  constitutive  model.  The  proposed  method  considered  the  rate  form  of  the 
constitutive  model  to  be  fitted  to  time-differentiated  stress-strain  data.  He  concluded  that 
there  is  a  definite  advantage  in  using  the  time-differentiated  stress-strain  data  to  fit  the 
model  since  this  form  of  data  tends  to  reduce  the  correlation  between  the  parameters  that 
define  the  model. 

Zhou  et  al.  (1998)  proposed  a  model  to  account  for  the  softening  effect  on  the 
yield  stress  produced  by  dynamic  recovery  and  heat  generation.  Their  viscoplastic  model 
consisted  of  a  hyperbolic-sine  strain  hardening  law  with  saturation  and  a  softening  term 
that  accounted  for  the  thermally  activated  effects  of  dislocation  motion.  The  model  was 
validated  with  torsional  tests  under  a  range  of  temperatures 

1.1.2  Fracture  Models 

The  papers  discussed  in  this  section  deal  with  the  modeling  of  dynamic  fracture. 
Most  of  the  articles  deal  with  ductile  fracture  in  the  form  of  void  nucleation  and  growth, 
but  there  are  some  review  papers  that  deal  with  the  dynamic  fracture  problem  in  general, 
including  brittle  fracture.  The  first  research  papers  that  are  discussed  do  not  deal  with 
fracture  but  with  the  deformation  of  porous  solids.  The  reason  for  presenting  a 
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discussion  of  these  papers  in  this  section  is  that  they  form  the  basis  for  the  latter 
development  of  void  growth  models  in  fracture  micro-mechanics. 

Mackenzie  (1949)  presented  a  model  to  calculate  the  effective  bulk  and  shear 
modulus  of  a  solid  containing  spherical  holes,  i.e.  material  micro-damage.  The  effective 
elastic  moduli  were  calculated  by  applying  a  stress  to  the  outer  spherical  boundary  of  a 
sphere  with  a  spherical  void  and  comparing  the  strains  produced  in  the  actual  solid  with 
the  strains  produced  by  the  same  stress  in  a  homogeneous  isotropic  voidless  solid. 

McClintock  (1968)  presented  a  simple  fracture  criterion  for  metals  driven  by  the 
growth  of  existing  holes  in  the  material.  His  model  was  based  on  the  growth  of 
cylindrical  holes  of  elliptical  cross  section  with  axes  parallel  to  the  principal  directions  of 
the  applied  stress. 

Herrmann  (1969)  proposed  a  constitutive  model  to  describe  the  compaction  of  a 
porous  solid  under  shock  loading.  The  theory  behind  the  model  was  phenomenological 
and  the  functions  that  define  the  model  had  to  be  fitted  to  empirical  data.  He  made  a 
series  of  simplifications,  including  the  assumption  that  there  is  no  plastic  shear 
deformation  and  that  the  compaction  process  is  insensitive  to  changes  in  temperature 
deviations  from  the  Hugoniot.  The  model  compared  well  with  experimental  data 
obtained  from  plate  impact  experiments. 

Carroll  et  al.  (1972)  derived  pore  collapse  relations  by  analysis  of  a  hollow  sphere 
under  an  external  pressure,  with  a  pressure-volume  relation  of  the  porous  material  similar 
to  the  one  proposed  by  Herrmann  (1969).  They  did  a  study  of  the  initial  elastic- void- 
volume  change,  the  transition  to  an  elastic-plastic  phase  and  the  final  plastic  phase.  They 
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concluded  that  the  first  two  phases  do  not  have  a  big  effect  on  the  compressibility  of  the 
porous  material.  The  plasticity  model  was  simplified  to  an  elastic-perfectly-plastic 
material  without  rate  effects,  and  thermal  effects  were  neglected. 

Butcher  et  al.  (1974)  presented  some  improvements  to  the  original  model 
proposed  by  Carrol  et  al.  (1972).  The  improvements  included  the  addition  of  work 
hardening  and  rate  dependence  through  a  viscosity  term.  The  predictions  of  this 
improved  model  were  compared  with  experimental  data  on  2024  Aluminum  with  good 
agreement. 

Shima  et  al.  (1976)  proposed  a  yield  criterion  for  porous  materials  that  coupled 
the  onset  of  yield  of  the  matrix  material  due  to  deviatoric  stresses  and  the  onset  of  yield 
on  the  pore  wall. 

Seaman  et  al.  (1976)  developed  models  of  ductile  and  brittle  fracture.  The 
models  included  the  nucleation  and  growth  of  voids  or  cracks  and  the  coalescence  of 
cracks.  The  nucleation  and  growth  models  were  phenomenological  and  included  a 
threshold  stress  and  material  viscosity  in  the  growth  model.  They  also  proposed  a 
modification  to  the  pressure-volume  relation  based  in  the  MacKenzie  (1949)  model  for 
the  bulk  modulus,  and  a  dependence  of  the  yield  stress  on  the  current  void  volume  and 
density. 

Gurson  (1977)  proposed  a  yield  function  for  porous  materials  and  a  void 
nucleation  model.  The  yield  function  was  used  as  a  plastic  potential  to  define  a  flow  rule 
that  is  independent  of  rate  and  thermal  effects,  but  considers  the  void  volume  fraction  in 
the  damage  material.  The  nucleation  model  is  based  on  the  statistical  distribution  of  solid 
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particles  in  the  matrix.  Thus,  a  void  is  nucleated  when  a  critical  mean  stress  is  reached 
and  the  particle  surface  is  debonded  from  the  matrix  material.  In  his  model  he  considered 
spherical  and  cylindrical  geometries  to  describe  the  shape  of  the  voids. 

Norris  et  al.  (1973)  proposed  a  phenomenological  model  that  describes  the  onset 
of  fracture  as  a  function  of  plastic-strain  and  mean-stress.  Fracture  occurs  when  a  critical 
value  of  a  cumulative-damage  parameter  is  reached  over  a  length  that  is  characteristic  of 
the  material.  They  tested  their  model  using  simulations  of  tests  that  included  simple 
tension,  circumferentially  notched  tension,  Charpy  V-notched,  and  pre-cracked  compact 
tension. 

Wilkins  et  al.  (1980)  also  proposed  a  cumulative-strain-damage  model  for  ductile 
fracture.  The  damage  parameter  was  defined  to  be  a  function  of  the  equivalent  plastic 
strain,  the  mean  stress  and  the  stress  deviator.  In  this  model  fracture  begins  when  the 
cumulative  damage  parameter  exceeds  a  critical  damage  value  over  a  critical  distance. 
The  model  was  calibrated  and  tested  with  simple  tensile,  cylinder  impact,  torsional  and 
notched  bar  tensile  tests, 

Johnson  (1981)  presented  the  complete  derivation  of  the  void  collapse  model 
presented  by  Carroll  et  al.  (1972)  with  significant  additions.  The  first  improvement 
proposed  by  Johnson  was  the  generalization  of  the  model  to  include  not  only  void 
collapse  but  also  void  distension.  He  expanded  on  the  work  by  MacKenzie  (1949)  to 
describe  the  elastic  moduli  dependence  on  void  volume  fraction  by  modifying  it  to  make 
the  shear  modulus  go  to  zero  as  the  void  volume  fraction  approached  unity.  In  addition, 
he  included  rate  effects  in  the  yield  function  through  a  viscosity  term,  but  left  out  thermal 
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softening  and  strain  hardening  effects.  Even  though  he  did  not  include  thermal  softening 
effects,  he  inferred  that  they  were  important  since  he  had  to  reduce  the  viscosity  of  the 
material  to  10P  from  room  temperature  values  of  approximately  ~105P  to  obtain  results 
comparable  to  experiments.  The  inertial  ef£ct  of  the  pore  wall  motion  was  neglected. 
He  also  addressed  the  problem  of  fragmentation  of  an  expanding  ring  through  the 
inclusion  of  a  heterogeneous  initial-void-volume-fraction  distribution  across  the  ring. 

Tvergaard  et  al.  (1982,  1984, 1986  and  1987)  and  Becker  et  al.  (1988)  presented  a 
constitutive  model  to  describe  ductile  fracture  in  polycrystalline  metals  by  nucleation, 
growth  and  coalescence  of  voids.  They  considered  a  material  with  voids  of  different 
orders  of  magnitude  in  which  the  material  has  an  initial  array  of  large  cylindrical  voids 
and  the  material  between  these  voids  is  modeled  through  a  constitutive  model  for  porous 
materials.  Thus,  fracture  occurs  when  the  large  voids  are  linked  by  micro-cracks,  which 
are  a  result  of  small-scale  nucleation,  growth  and  coalescence  of  micro-voids.  Tvergaard 
(1982)  proposed  a  yield  condition  for  a  porous  solid  that  was  based  on  the  work  of 
Gurson  (1977),  with  modifications  to  improve  the  modeling  of  shear  band  localization. 
Tvergaard  and  Needleman  (1984)  also  proposed  a  modification  to  the  yield  condition  to 
account  for  failure  through  the  introduction  of  a  critical  and  ultimate  void-volume- 
fraction.  They  also  proposed  phenomenological  nucleation  and  growth  models  that  are 
strain  and  stress  driven.  Tvergaard  (1996)  used  his  model  to  study  the  effect  of  void  size 
distribution  in  the  matrix  material  on  the  fracture  process.  He  concluded  that  for  large 
void  volume  fraction  the  small  voids  grow  at  a  higher  rate  and  conversely,  for  small  void 
volume  fraction  the  large  voids  grow  at  a  higher  rate. 
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Hancock  et  al.  (1985)  discussed  the  results  of  a  finite  element  solution  for  a 
spherical  elastic  inclusion  in  a  plastically  deforming  matrix.  These  results  were  used  in 
conjunction  with  some  experimental  work  using  notched  tensile  specimens  to  generate 
multiaxial  states  of  stress  from  which  the  local  conditions  leading  to  decohesion  (void 
nucleation)  of  the  inclusion/matrix  interface  were  determined.  One  important  feature  of 
this  study  was  the  conclusion  that  there  is  a  statistical  distribution  of  interfacial  strength. 
Some  of  the  voids  had  a  growth-threshold  mean-stress  approximately  seven  times  the 
yield  stress  and  some  weakly  bonded  particles  behave  as  pre-existing  voids. 

Carroll  (1985)  and  Kim  and  Carroll  (1987)  presented  an  extension  to  his  work  of 
(Carroll  1 972)  to  include  a  study  of  different  strain  hardening  laws  that  allowed  for  a 
closed-form  solution  to  the  variable  mean  stress  threshold  for  void  volume  fraction 
growth.  The  strain  hardening  laws  included  an  elastic  response,  two  saturation  hardening 
laws,  a  pseudo-exponential  hardening  and  a  power-exponential  hardening  law.  He  also 
included  some  discussion  on  the  behavior  of  these  response  laws  in  terms  of  the  shape  of 
the  threshold  mean  stress  as  a  function  of  void  volume  fraction.  In  addition,  Carroll  et  al. 
(1986)  presented  a  study  of  the  effect  of  temperature  on  the  viscoplastic  response  of  a 
single  void  under  compressive  mean  stress,  concluding  that  the  viscoplastic  response  is 
better  modeled  when  the  temperature  effect  is  included.  This  was  done  by  introducing  a 
temperature  dependent  yield  stress  and  viscosity  model,  and  an  energy  rate  balance 
equation. 

Perzyna  (1986)  presented  a  model  for  dynamic  fracture  that  was  based  in  part  on 
the  work  of  Carroll  (1972,  1985)  and  Johnson  (1981).  Perzyna  simplified  the  void 
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growth  model  by  also  excluding  the  inertial  effects,  which  transformed  the  second  order 
nonlinear  differential  equation  for  void  growth  to  a  first  order  differential  equation.  He 
also  proposed  a  viscoplastic  model  with  linear  work  hardening  for  the  void  growth  model 
but  left  out  all  temperature  effects.  The  most  important  contribution  by  Perzyna  was  the 
formulation  of  an  elastic-viscoplastic  model  of  the  voided  solid  that  inherently  described 
the  dilatational  effects  and  fracture  generated  by  the  nucleation,  growth  and  coalescence 
of  microvoids.  His  formulation  was  simpler  than  the  one  proposed  by  Tvergaard  and 
Needleman  (1984). 

Curran  et  al.  (1987)  presented  a  thorough  review  of  proposed  models  that 
describe  the  dynamic  failure  mechanisms  in  polycrystalline  solids  from  micro  and 
macroscopic  points  of  view.  In  addition,  they  also  presented  a  brief  discussion  of  the 
experimental  techniques  used  to  characterize  the  dynamic  fracture  process.  They 
presented  a  discussion  of  the  models  for  void  and  crack  nucleation,  a  discussion  of  the 
different  models  used  to  describe  the  growth  of  microvoids  or  shear  bands  and  their 
interaction  leading  to  the  coalescence  of  microvoids  or  cracks. 

Grady  (1988)  discussed  the  spall  strength  of  solids  in  terms  of  its  fracture  energy. 
His  study  had  the  objective  of  identifying  which  properties  of  the  material  influence  the 
damage  process  in  both,  brittle  and  ductile  failure.  He  defined  the  spall  strength  for 
brittle  fracture  to  be  controlled  by  the  fracture  toughness,  and  for  ductile  fracture  by  the 
yield  stress.  Due  to  the  simplifications  and  assumptions  in  the  modeling,  he  reached 
conclusions  that  have  been  proven  to  be  inaccurate,  such  as  the  conclusion  that  ductile 
spall  strength  is  rate  independent. 
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Rajendran  et  al.  (1989)  presented  a  dynamic  failure  model  based  on 
phenomenological  observations.  They  modeled  the  plastic  response  of  the  matrix 
material  with  the  Bodner  and  Partom  model,  which  does  not  consider  a  yield  condition, 
and  the  aggregate  with  voids  with  a  model  similar  in  form  to  Gurson’s  (1977)  yield 
model.  They  tested  their  model  with  plate  impact  experiments  of  OFHC  Copper  at 
relatively  low  velocities  (<200m/s). 

In  a  series  of  papers  from  1989  to  1996  Eftis  and  Nemes  generalized  Perzyna’s 
viscoplastic-damage  constitutive  model  and  calculation  of  the  microvoid  growth  rate,  by 
incorporating  nonlinear  strain  hardening  that  saturates  for  the  viscoplastic  material 
behavior.  This  elastic-viscoplastic  microdamage  constitutive  model,  appropriate  for 
isothermal  ductile  fracture  of  polycrystalline  materials,  was  used  to  successfully  describe 
the  main  features  of  plate  impact-induced  spall  fracture  for:  (i)  Rectangular  plates,  (ii) 
oblique  impact  of  rectangular  plates,  (iii)  circular  plates  with  the  flyer  having  smaller 
diameter  then  the  target,  and  (iv)  very  high  strain  rate  fracture  of  unnotched  round  tensile 
bars.  Th£  ionstitutive-microdamage  equations  have  also  been  employed  to  analytically 
demonstrate  the  known  experimental  fact  that  elasto-viscoplastic  waves  can  be 
propagated  through  damaged  softened  (i.e.  unstable)  polycrystalline  materials. 

Cortes  (1992)  presented  a  pore  growth  model  that  was  based  on  the  model 
proposed  by  Carroll  (1972,  1985)  and  Johnson  (1981)  with  the  inclusion  of  linear 
hardening  and  viscous  and  thermal  effects  but  with  a  simplified  definition  of  strain.  He 
also  proposed  a  modification  to  the  original  model  to  include  the  effects  of  deviatoric 
stresses  on  the  growth  of  microvoids.  He  concluded  that  the  thermal  effect  on  the  void 
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volume  fraction  evolution  was  negligible,  which  contradicts  the  findings  of  Johnson 
(1981)  and  Carroll  and  Kim  (1986).  Cortes  et  al.  (1995)  modified  his  model  to  include 
void  nucleation  and  modified  the  growth  threshold-mean-stress  using  Gurson’s  (1977) 
yield  model.  He  used  his  new  model  to  simulate  spall  fracture  in  plate  impacts  and 
explosive  loading  of  plates  with  good  agreement  between  experiment  and  simulation 
results  in  terms  of  the  velocity  profile  of  the  target  rear  surface. 

Benson  (1993)  presented  a  study  of  the  effect  of  void  distribution  on  the  dynamic 
growth  and  coalescence  of  voids.  Computer  simulations  were  made  of  different  patterns 
of  cylindrical  voids  in  OFHC  Copper  and  4340  Steel.  The  Johnson-Cook  model  was 
used  to  describe  the  viscoplastic  behavior  of  the  material  and  a  simple  limit-strain 
fracture  model  was  implemented.  The  most  important  conclusion  from  his  work  was  that 
the  fracture  path' for  the  two  materials  was  different.  While  for  OFHC  Copper  fracture 
was  reached  when  the  voids  grew  until  they  connected,  4340  fractured  by  failure  of  the 
material  between  the  voids.  Benson  concluded  that  the  difference  in  fracture  paths  was 
due  to  the  large  difference  in  the  strain-threshold  for  fracture. 

Tong  et  al.  (1993)  presented  a  study  of  the  inertial,  thermal  and  viscous  effects  on 
the  void  collapse  obtained  through  a  model  similar  to  the  Carroll  and  Johnson  models  and 
presented  comparative  results  with  Carroll’s  model.  In  1995  Tong  did  a  similar  study  of 
the  inertial  effects  on  the  void  collapse  predicted  by  the  Gurson  yield  model.  In  all  cases, 
he  concluded  that  inertial  and  dynamic  effects  in  general  have  a  significant  effect  on  the 
pore  collapse  predictions. 

Addessio  et  al.  (1993),  Lee  et  al.  (1994,  1996)  presented  studies  of  some  of  the 


models  already  mentioned  in  this  review  for  specific  applications  (i.e.  compaction  of 
cylinders,  formability,  etc)  and  proposed  small  improvements  with  suggested 
shortcomings  and  strengths  of  the  models. 

Kanel  et  al.  (1996)  presented  the  results  of  a  series  of  experiments  done  on 
Aluminum  ADI  and  Magnesium  Mg95  plates  to  investigate  the  effect  of  temperature  on 
the  spall  strength  of  these  materials.  The  conclusion  was  that  a  higher  initial  temperature 
of  the  specimens  caused  the  spall  strength  to  decrease,  which  again,  contradicts  the 
conslusion  of  Cortes  and  supports  the  conclusions  by  Carroll  and  Kim  (1986)  and 
Johnson  (1981). 

1.1.3  Computer  Codes 

The  literature  review  in  terms  of  the  numerical  algorithms  and  computer  codes 
was  not  as  extensive  as  in  the  areas  of  constitutive  models  and  fracture,  and  in  the 
following  section  on  materials  characterization  of  after  shock  specimens.  The  objective 
of  the  review  on  computer  codes  is  directed  towards  providing  an  overview  of  the 
numerical  algorithms  used  in  the  commercial  and  public  domain  codes  used  to  simulate 
hypervelocity  impacts,  and  not  to  give  a  detailed  exposition  of  the  evolution  of  every 
component  and  algorithm  that  compose  current  computer  codes.  The  reason  for  this  is 
that  the  research  work  emphasized  the  development  of  an  elastic-viscoplastic  micro- 
damage  constitutive  model  appropriate  for  hypervelocity  impact  and  its  implementation 
as  separate  subroutines  in  an  existing  computer  code. 

Current  computer  codes  and  their  algorithms  are  based  on  the  work  by  Wilkins  et 
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al  (1974)  where  they  described  the  details  for  the  development  of  a  three-dimensional 
finite  element/finite  difference  computer  program  to  solve  high-shock  and  high  non¬ 
linear  deformation  problems  in  a  Lagrangian  grid. 

Liu  et  al.  (1986)  and  Benson  (1989)  presented  the  conservation  laws,  the 
constitutive  relations,  and  the  pressure-volume  relations  for  the  modeling  of  path- 
dependent  materials  in  an  arbitrary  Lagrangian-Eulerian  finite  element  method. 

Johnson  et  al.  (1987)  introduced  an  eroding  algorithm  to  be  incorporated  in  a 
Lagrangian  finite  element  code  to  simulate  projectile  penetration/perforation  of  thick 
plates.  The  objective  of  this  eroding  algorithm  is  to  eliminate  computational  cells  when 
they  reach  a  limit  strain. 

McGlaun  et  al.  (1990)  gave  a  brief  description  of  the  models  and  features 
available  in  the  CTH  Eularian  code.  Trucano  et  al.  (1990)  showed  the  results  of  CTH 
simulations  of  hypervelocity  impacts.  They  reported  on  the  effect  of  grid  density, 
advection  and  interface  tracking  algorithms  and  other  factors  on  the  simulation  errors. 

Anderson  (1987)  provides  an  extensive  review  of  the  evolution  of  “Hydrocodes”. 
He  discusses  Lagrangian  and  Eularian  formulations,  including  advantages  and 
disadvantages,  the  concept  of  artificial  viscosity  to  deal  with  discontinuities  of  shocks, 
and  finally  the  treatment  of  constitutive  models,  pressure-volume  relations,  eroding 
algorithms  and  fracture. 

1.1.4  Microstructural  Studies 

This  section  of  the  review  will  briefly  discuss  microstructural  studies  as  they 
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relate  to  hypervelocity  impact  specimens.  The  papers  discussed  in  this  section 
concentrate  on  the  effects  of  pre-impact  microstructure  and  the  microstructural  features 
observed  in  after-impact  specimens. 

Meyers  et  al.  (1981)  presented  a  short  review  of  the  mechanisms  that  lead  to 
defect  generation  under  shock  deformation.  They  discussed  the  formation  of  point,  line 
and  planar  defects,  and  phase  transformations. 

Gray  III  (1993)  gives  a  review  of  the  effects  of  shock  deformation  on  material 
response.  He  discusses  the  techniques  and  design  parameters  used  in  shock  recovery 
experiments,  and  the  influence  of  shock  parameters,  including  peak  pressure,  pulse 
duration  and  loading  rate  on  post-shock  structure  or  properties  of  the  material. 

Ferreyra  el  al.  (1995)  provided  a  study  of  the  effects  of  target  microstructure  on 
crater  depth,  diameter  and  depth  to  diameter  ratios  in  a  Copper  plate  impacted  by  a 
spherical  projectile  at  high  velocity.  They  conclude  that  there  is  a  profound  effect  of 
target  microstructure  on  cratering  characteristics. 

Rivas  el  al.  (1995)  gave  a  brief  discussion  of  the  different  after-shock 
microstructural  characteristics  observed  in  OFHC  Copper,  1100  and  6061-T6  Aluminum 
targets. 

Murr  at  al.  (1996)  discussed  the  deformation  processes  and  microstructures 
associated  with  different  high  deformation  configurations.  They  discussed  the 
microstructural  features  and  deformation  mechanisms  associated  with  plate  impact 
experiments,  shaped  charges,  explosively  formed  penetrators,  spherical  projectile 
hypervelocity  impacts  and  long  rod  penetration  of  thick  targets. 
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Quinones  (1996)  provided  a  very  thorough  study  of  the  microstuctural  features 
observed  in  post  impact  OFHC  Copper  plates.  The  study  included  the  discussion  of 
microstructural-features  and  hardness  profiles  across  the  thickness  of  the  target  plate  after 
impact  at  different  velocities.  She  also  simulated  the  imparls  on  a  computer  code  and 
correlated  the  hardness  profile  with  residual  yield  stress  as  defined  by  the  Johnson-Cook 
yield  model.  Quinones  et  al.  (1998)  extended  the  simulation  study  to  include  the 
simulation  of  spall  fracture  using  a  limit  stress  fracture  model. 

Murr  et  al.  (1997)  and  Ferreyra  et  al.  (1997)  provided  a  study  of  microstructure 
after  high  velocity  and  hypervelocity  impact  cratering  of  Copper  targets.  They  conclude 
that  the  target  microstructure  has  a  strong  effect  on  the  geometrical  characteristics  of  the 
craters,  especially  the  dislocation  structure  while  grain  size  seems  to  have  an  insignificant 
effect. 


1.2  Objective 

The  research  presented  in  this  dissertation  seeks  to  expand  the  capability  of 
current  constitutive  models  for  the  simulation  of  hypervelocity  impact  problems.  This  is 
accomplished  through  the  development  of  a  more  general  constitutive-microdamage 
model  that  can  describe  high  shock  compression,  dilatation,  shear  and  material 
microdamage  evolution  leading  to  fracture,  assuming  thermoelastic-thermoviscoplastic 
material  behavior.  The  microdamage  model  presented  here  incorporates  a  non-linear 
pressure  volume  relation  for  shock  compression  and  dilatation,  along  with  expressions  for 
the  elastic  moduli  that  include  volume  strain  dependence  as  well  as  thermal  and 
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microdamage  softening.  The  viscoplastic  response  equations  incorporate  nonlinear  strain 
and  strain  rate  hardening  as  well  as  thermal  and  microdamage  softening  with  temperature 
dependent  viscosity,  and  include  plastic  compressibility  due  to  microvoid  volume 
fraction  change.  The  melt  temperature  is  assumed  to  be  volume  strain  dependant.  The 
microdamage  evolution  model  includes  microvoid  void  nucleation  with  a  variable  mean 
stress  threshold  and  microvoid  growth  with  void  interaction  effects.  The  microvoid 
growth  model  incorporates  a  viscoplastic  material  model  with  nonlinear  strain  and  strain 
rate  hardening  as  well  as  thermal  softening,  and  a  temperature  dependant  viscosity. 

The  work  presented  in  this  dissertation  improves  upon  previous  works,  related  to 
hypervelocity  impact  modeling  and  simulation,  in  several  key  issues. 

1)  The  constitutive  model  presented  here  is  the  most  extensive  in  the  sense  that  it 
incorporates  more  of  the  known  polycrystalline  material  behavior  issues,  which 
were  described  in  this  section. 

2)  In  the  constitutive  model  presented  here,  the  softening  effect  of  microdamage  is 
incorporated  into  the  Johnson-Cook  model  for  the  yield  function,  as  wed  as  in  the 
equations  for  the  elastic  moduli. 

3)  In  the  microvoid  nucleation  rate  equation,  the  threshold  mean-stress  is  not 
assumed  constant  but  is  expressed  as  a  function  of  equivalent  plastic  strain  and 
strain  rate. 

4)  A  temperature  rate  equation  local  to  the  wall  of  a  growing  microvoid  is 
formulated.  The  temperature  calculated  with  this  equation  is  substituted  in  the 
expressions  for  the  temperature  dependent  material  viscosity  and  threshold  mean- 
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stress  for  microvoid  growth.  This  provides  analytical  justification  for  the  reduced 
temperature-independent  viscosity  used  by  Johnson  (1981)  and  Eftis  et  al.  (1991), 
and  represents  an  important  improvement  over  previous  models. 

5)  In  the  expression  for  the  viscoplastic  rate  of  deformation  tensor,  the  volumetric 
component  (due  to  microvoid  volume  change)  is  not  explicitly  coupled  to  the 
deviatoric  component.  This  corrects  the  problems  associated  with  the  Perzyna’s 
model  (Perzyna,  1986)  that  give  a  non-zero  value  to  the  volumetric  contribution 
from  micro  void  growth  even  when  the  micro  void  volume  rate  of  change  is  zero. 

6)  Most  of  the  constitutive  models  for  the  simulation  of  hypervelocity  impact 
problems  presented  in  the  literature  are  tested  with  only  one  computer  simulation, 
while  the  simulations  reported  in  this  dissertation  (see  Chapter  6)  cover  a  wide 
range  of  plate  target  thickness  that  result  in  widely  varied  levels  of  damage. 
These  simulations  provide  a  measure  of  the  capability  and  performance  of  the 
microdamage  constitutive  model  in  the  simulation  of  hypervelocity  impact 
problems  under  varied  geometrical  configurations. 


CHAPTER  2 


f 

•  RANKINE-HUGONIOT  RELATIONS 


Extensive  experimental  characterizations  of  the  pressure-volume  (P-V)  relations 
for  various  materials  for  a  wide  range  of  shock  pressures  have  been  developed  since  the 
1950’s  (Walsh  et  al.  1957).  Most  of  the  data  compiled  has  been  produced  through 
contact  explosive  detonation,  or  the  “flyer  plate”  experiment  where  a  projectile  plate  is 
explosively  accelerated  towards  a  target  plate  in  a  configuration  that  produces  planar 
shock  waves  across  the  projectile  and  target  plates.  The  data  is  usually  presented  as  P-V 
curves  or  curves  that  relate  the  shock  wave  speed  Us  to  the  particle  velocity  UP  behind 
the  shock  front.  It  has  been  discovered  experimentally  that  for  polycrystalline  metals 
under  a  wide  range  of  shock  pressures  the  Us-UP  curve  can  be  expressed  as  a  linear 
equation, 

Us=c0  +  sUp  (2.1) 

where  c0  is  the  sound  speed  in  the  material  at  atmospheric  pressure  and  room  temperature 
(~300°K) ,  s  is  the  slope  of  the  US-UP  linear  relation  and  UP  is  the  particle  speed  behind 
the  shock  front.  It  is  important  to  note  that  the  experimental  setup  produces  planar  shock 
waves  where  the  deviatoric  stresses  are  non-zero  and  thus  some  plastic  deformation  and 
“residual”  temperature  is  to  be  expected  if  the  dynamic  yield  stress  is  exceeded.  Figure 
2.1  shows  experimental  data  points  for  Copper  in  the  Us-UP  plane  with  a  linear  fit 
through  the  points. 
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Figure  2.1  High-pressure  Hugoniot  data  for  Copper  obtained  by  several 
investigators.  The  numerals  along  the  curve  indicate  pressures  in  GPa.  The  line  on  the 
small  inset  represents  the  entire  range  of  the  larger  figure,  while  the  three  points  represent 
the  results  of  Soviet  ultrahigh  pressure  experiments  (Davison  et  al.  1979). 

The  development  of  a  P-V  equation  for  shock  loading  requires  the 
implementation  of  the  “jump  conditions”  across  the  shock  wave  coupled  with  the  Us-Up 
linear  relation.  The  “jump  conditions”  are  determined  by  the  conservation  of  mass,  linear 
momentum  and  energy  that  relate  the  density,  particle  velocity  and  stresses  across  the 
discontinuity  produced  by  the  moving  shock  wave,  and  can  be  stated  at  first  in  general 
terms  by  the  well  known  jump  relations  (Figure  2.2), 

P*  (un  ~  vn )  =  P"  (u„  -  v" ) ,  conservation  of  mass  (2.2) 

p*v*(un-v;)+r.n 

/  \  ,  conservation  of  momentum  (2.3) 

=  p-v'(un-v-)+T--n 
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^Iv+.v++E+j(un-v;)+t+-v+-q-n 
f\v  ■ v'  +  E"](un  -  V')+ 1'  •  V-  -q  •  n 


,  conservation  of  energy 


(2.4) 


where  p  is  the  density,  un  is  the  speed  of  the  shock  wave  along  a  unit  vector  normal  to  the 
shock  front  (n),  vn  is  the  particle  speed  along  vector  n,  v  is  the  particle  velocity  vector,  T 
is  the  Cauchy  stress  tensor,  t  is  equal  to  Tn,  which  is  the  traction  vector  across  the  wave 
front,  E  is  the  internal  energy  per  unit  mass  and  q  is  the  heat  flux.  In  these  equations,  the 
functions  (<f>+)  and  (<)>')  indicate  the  limits  of  the  function  4>  as  the  shock  waye  is 
approached  from  the  front  (material  into  which  the  wave  is  advancing)  and  from  the 
back,  respectively. 


u 
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Figure  2.2  Depiction  of  the  state  of  the  material  ahead  and  behind  the  Shock-Front. 


If  we  make  the  shock  front  a  planar  surface  moving  in  the  xi  direction  of  a 
rectangular  Cartesian  coordinate  system  so  that  |n|  =  n,  =1,  q=0  (adiabatic  conditions), 
and  denote  un=ui  then  Equations  (2.2)  -  (2.4)  can  be  expressed  as, 
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P+(u- v;)= p-(u- vr)  (2.5) 

px(u  -  vr )+  V\  =  PX(U  -  V- )+  Tn  (2.6a) 

T,+2  =  T,-  (2.6b) 

TI3  =  T13  (2.6c) 

P'^ky+E-jtu-vrJ+T'v^p-f^v-y  +  E-Vv^+T^vr  (2.7) 


Equations  (2.6b)  and  (2.6c)  state  that  the  shear  stress  fields  T12  and  T13  are  continuous 
across  the  planar  shock  wave.  If  we  now  let  v+=0  and  make  u=Us,  vj-  =  Up,  p+=p°, 

p'=p,  making  use  of  Equations  (2.1)  and  (2.5)  -  (2.7)  we  can  simplify  the  jump  equations, 
and  obtain  the  Rankine-Hugoniot  equations 


P0US  =  P(us  -  Up) 

(2.8) 

T|i  —  T|i  =  P0UPUS 

(2.9) 

.k+Tjfi-e/) 

E"  -  E*  -  1  V  /Po) 

2  P 

(2.10) 

Combining  Equations  (2.1),  (2.1 1)  and  (2.12)  yields, 

nr  +  y-  PoC0  V 

"  "'[1-srf 

(2.11) 

where,  vj/  =  1  -  — ,  is  the  volume  strain. 

n 

(2.12) 

Equation  (2.11)  represents  the  change  in  stress  Tn  from  an  initial  state  to  a 


shocked  state.  This  equation  can  be  modified  and  generalized  to  a  shock  wave  front  with 
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arbitrary  orientation.  Then  the  uniaxial  stress  can  be  replaced  by  the  negative  of  the 
mean  stress  (denoted  as  a  ‘pressure’)  P  =  trT  so  that  Equation  (2.1 1)  can  be  stated  as, 


p  _  p  —  PocoV  —  p 
H  0  [l-syf  H‘ 


(2.13) 


where  Ph  is  the  Hugoniot  pressure  above  a  zero  pressure  reference,  Po  is  the  pressure 
acting  on  the  material  ahead  of  the  shock  front  above  the  zero  pressure  reference  and 
PHu  is  the  Hugoniot  pressure  above  Po.  When  the  volume  strain  \j/  is  positive  signifying 

volume  reduction,  the  ‘pressure’  P  is  positive,  and  conversely  when  the  volume  strain  is 
negative  signifying  dilatation,  P  is  negative. 

Equation  (2.13)  is  generally  called  the  Hugoniot  pressure  and  represents  the  locus 
of  all  the  shocked  states  as  depicted  in  Figure  2.3.  When  a  material  is  shocked,  the 
pressure  does  not  follow  the  Hugoniot  (P-vj/)  pressure,  but  increases  discontinually  from 
its  initial  state  Po  to  Pi  following  the  Rayleigh  line  (Davison  et  al.  1979  and  Boslough  et 
al.  1993).  Similarly,  from  Equation  (2.10)  the  change  in  internal  energy  due  to  the  shock 
compression  can  be  shown  to  be, 


_  _  1  (p  +  Pok  _p 

lH  l0",  “ 

2  Po 


(2.14) 


where  Eh  is  the  Hugonion  energy  above  a  zero  energy  reference,  Eo  is  the  energy  of  the 
material  ahead  of  the  shock  front  above  the  zero  energy  reference  and  EHo  is  the 

hugoniot  energy  above  Eo.  It  should  be  noted  that  Equation  (2.14)  gives  the  area  between 
the  Rayleigh  line  and  the  zero  reference  pressure  (Figure  2.3). 
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Figure  2.3  Shock  Hugoniot  pressure  curve  and  Rayleigh  line  for  a  given  pressure  Pi. 


Equation  (2.13)  and  (2.14)  represent  the  pressure  and  energy  increments  associated  with  a 
change  in  volume  strain  (v|/),  but  do  not  include  thermal  contributions  from  plastic  work 
associated  with  the  deviatoric  stresses. 


CHAPTER  3 


GRUNEISEN  EQUATION  OF  STATE 


The  thermal  pressure  within  a  material  is  associated  with  its  thermal  energy, 
which  can  be  calculated  approximately  by  assuming  the  ion  cores  of  a  crystal  lattice  as 
N-many  simple  harmonic  oscillators  (the  oscillations  of  the  free  electron  gas  are  not 
considered).  The  oscillators  are  vibrating  in  the  normal  modes  of  the  three-dimensional 
system  at  frequencies  Vj(u)  that  are  assumed  to  depend  only  upon  the  specific  volume 
(Walsh  et  all.  1957).  From  quantum  mechanical  calculations  of  the  energy  of  a  system  of 
3N  oscillators  coupled  with  a  statistical  mechanics  representation  of  the  thermal  energy 
of  the  lattice,  for  any  temperature  T>0°K  the  total  energy  can  be  shown  to  be  given  by 


E(»,T)-*0(.)+i!>,+f;  ^ 


ftv./ 

i  e  At_1 


(3.1) 


.  .  1  3N 

where,  Q(u)  is  the  binding  energy  of  the  lattice  at  0°K,  -^/zv,  is  the  ground  state 


3N  hv 

energy  of  the  0°K  oscillations  at  the  lowest  energy  state  and  ]>]-— i —  ls  the  thermal 

i  Q  /k T  _  | 

energy  of  the  oscillating  lattice,  with  the  summation  being  over  the  3N  normal  modes  of 
the  crystal.  The  constants  h  and  k  are  the  Planck  and  Boltzman  constants  respectively. 
The  Helmholtz  free  energy  can  be  expressed  as 

H  =  E  -  TS  .  (3.2) 

S  is  the  thermodynamic  entropy  defined  from  statistical  mechanics  as 
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S  =  KlnYeE(n)/‘T+- 


where 


E(n)  =  O0+^(n  +  j/J^vi  . 


and  n  is  the  occupation  number  of  the  phonon  energy  levels  oscillating  at  the  r 
frequencies.  Substituting  Equation  (3.3)  in  (3.2)  leads  to, 


H  =  -tcTln 


E(n)/*T 


Then,  substituting  Equation  (3.4)  in  (3.5)  leads  to  the  following  expression  for  the 
Helmholtz  free  energy 


1  3n  jn  r  M, 

H  =®0  +— +KT]T]ln  1-e  'v 

2  i  i  - 


The  external  pressure  is  defined  as  the  specific  volume  derivative  of  the  free  energy  at 
constant  temperature  taken  as  negative  for  volume  compression,  giving  the  result 


dO  A 


\  d  u  JT  2  \  u  Y 


e  /kT-l 


where 


Y,(v)_  if  din 


1  I  dv, 


\)  dim)  Jr  v,  ^  dv  )T 


Equation  (3.8)  represents  what  is  known  as  the  Gruneisen  ratio  at  the  ith  frequency  and  is 
a  function  of  the  specific  volume  (u).  A  general  assumption  that  is  used  to  simplify 
Equation  (3.7)  is  to  take  all  the  frequencies  (vj)  for  the  3N  harmonic  oscillators  as  being 


equal.  This  leads  to  the  following  simplification  for  the  Gruneisen  parameter 
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YiW  =  YiW.  (3.9) 

Substituting  Equation  (3.9)  in  (3.7)  leads  to  the  following  simplified  definition  of  the 


pressure, 


P  =  - 


^dO0"|  ,  1  YMyfo.  YHy  H 

2  v  j  ‘  v 


v  du  yT 


(3.10) 


-1 


If  we  now  define 


and 


1  3N 

e(u)0.k  =  ®0W+-Z*v.  ’  ^  ccdd  ener§y 

2  i 


3N  ft\> 

J — ,  as  the  thermal  energy, 
■  e  At_i 


(3.11) 


(3.12) 


then 


^d<D  ^ 


v  du  , 


y(d) 


P(v)o° k  - - -  +  -^-^E0(u) ,  represents  the  cold  pressure. 


u 


(3.13) 


By  combining  Equations  (3.10),  (3.12)  and  (3.1)  the  pressure  at  specific  volume  (u)  and 
temperature  (T)  can  be  expressed  as 


P(v,T)=P(«)t.1<+I^Ejtl,T). 

1) 


(3.14) 


This  definition  of  the  pressure  can  be  traced,  in  Figure  3.1,  by  following  the  0°K  isotherm 
up  to  specific  volume  (u)  or  volume-strain  (vy)  at  pressure  P(vy,0°K)  and  then  adding,  at 
constant  volume,  the  thermal  contribution  up  to  P(\|/,T). 

From  Equations  (3.1),  (3.1 1)  and  (3.12)  we  can  define  the  thermal  energy  as  follows: 
Elh(u,T)=E(u,T)-E(u)0.IC 


(3.15) 
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Substituting  Equation  (3.15)  in  (3.14)  we  reach  the  following  result, 


P(u,T)=  P(»)„  +lfcl(E(>.,T)-E(«)t.K) 

V 

(3.16a) 

P(U’T)-P(u)0.k  =^(E(u,T)-E(u)0...). 

(3.16b) 

Equation  (3.16a)  is  known  as  the  general  Mie-Gruneisen  Equation  of  state  (EOS). 
Equation  (3.16b)  represents  the  change  in  pressure  above  the  0°K  isotherm  (P-V-E) 
curve.  This  form  of  the  Mie-Gruneisen  EOS  makes  it  easy  to  visualize  the  fact  that  it  can 
be  rewritten  in  terms  of  any  other  (P-V-E)  curve.  Rewriting  Equation  (3.16a)  in  terms  of 
the  Hugoniot  curve  taken  as  reference, 

P(«,T)=Ph+^0(E(U,T)-Eh).  (3.17) 

t) 

This  definition  of  the  pressure  can  be  traced,  in  Figure  3.1,  by  following  the  Hugoniot 
curve  from  the  Po  point  up  to  Ph  and  then  adding  the  partial  thermal  contribution  by  the 
difference  between  the  total  energy  E(vg,T)  and  the  Hugoniot  energy  Eh- 
Then,  substituting  Equations  (2.13)  and  (2.14)  in  Equation  (3.17)  and  rearranging  terms, 
P(v,T)-P0  =PHo  +y(v)p(EHJ+Y(\)/)p(E(u,T)-E0)  (3.18a) 

P(V-T)-  P,  =  J)  +  7Mp(E(V;T)-E0).  (3.18b) 

Equation  (3.18a)  represents  the  change  in  pressure  above  the  reference  point  (Po,Eo) 
which  has  to  be  the  room  conditions  of  300°K  temperature  and  atmospheric  pressure, 
since  these  are  the  conditions  under  which  the  Us-Up  curves  (Eq  2.1)  that  form  the  basis 
for  the  Hugoniot  curve,  are  experimentally  obtained.  Thus  the  Mie-Gruneisen  (P-V-E) 


equation  can  be  simplified  as  follows  with  the  understanding  that  the  pressures  and 
energies  represent  changes  with  respect  to  the  reference  point  previously  mentioned, 

P('i'.T)= y)) + tMp  e(v,t),  v>o  (3.19) 

where  E(vy,T)  is  the  internal  energy  per  unit  mass  that  can  be  calculated  by  the  energy 
rate  balance  equation  appearing  in  Section  4.5. 

Figure  3.1  presents  a  qualitative  illustration  of  all  the  different  terms  that  lead  to 
the  formulation  of  the  Mie-Gruneisen  (P-y)  curve. 

E 

e(v|/,T) 

Eh 

e(v,300°K.) 

E(vy,0°K) 

e°(m/o,300°K) 

E(yo,0°K) 

M'°  V(P.T)  V 

Figure  3.1  Mie-Gruneisen,  Hugoniot,  0°K  and  300°K  Isothermal  Pressure- 
Volume  Strain  Curves. 
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The  volume  strain  dependence  of  the  Gruneisen  coefficient  y(\y)  of  Equation 
(3.19)  can  be  expressed  as  a  polynomial  with  empirically  determined  coefficients  (Walsh 
et  al.  1957).  Then,  to  a  first  order  approximation  y(vpr)  is  expressed  as 


yM=y0  +a 


V 

1-Vy 


(3.20) 


where,  for  example  the  coefficients  y0  and  A  have  the  values  (2.04,  -3.296),  respectively 
for  Copper  and  (2.13,  -7.245)  for  Aluminum-24  ST. 

Insertion  of  (3.20)  into  (3.19)  gives  the  nonlinear  pressure-volume  strain  relation 
(equation  of  state)  associated  with  shock  compression  of  polycrystalline  materials. 


P(V,T)  = 


PocoV 


[l-svpf 
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(  > 

r  ...  Vi 
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f  M 

V 

+A 

> 

4- 

Yo+A 

V 

1 

l1-vj 

ll-vj 

. 

/ 

U-vJ_ 

pE(\j/,T),  \\)>0  (3.21) 


It  will  be  recalled  that  in  the  development  of  Equation  (2.19)  the  ‘pressure’  term  P 


actually  represents  the  negative  of  the  mean  stress  P  =  ~~trT  =  -am ,  as  designated 

above,  E  is  the  internal  energy  per  unit  mass,  p,  p0  are  the  mass  densities  in  the  current 
and  reference  states,  (co)  is  the  sound  speed  and  (s)  is  the  experimentally  determined 
slope  in  the  shock  wave  speed-particle  speed  relation  (2.1) 
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CHAPTER  4 

CONSTITUTIVE  -  MICRODAMAGE  MODEL 

4.1  Introduction 

The  governing  balance  Equations  of  continuum  mechanics,  expressed  in  local 
form,  state  that  for  any  deforming  material  body  the  balance  of  mass,  linear  momentum, 
moment  of  momentum  and  energy,  respectively,  requires  that  the  following  Equations  be 


satisfied  at  each  point  of  the  body 

+  grad  p  •  v  +  p  div  v  =  0 

(4.1) 

dv 

divT  +  pb  =  p — 
dt 

(4.2) 

T  =  Tt 

(4.3) 

de 

p—  =  tr(T  -D)  +  ph-divq 

(4.4) 

where  p  is  the  mass  density,  v  is  the  velocity  vector,  T  the  Cauchy  stress  tensor  where  the 
superscript  denotes  transpose,  b  is  the  body  force  per  unit  mass,  D  is  the  rate  of 
deformation  tensor,  h  is  the  heat  source  per  unit  mass  and  q  is  the  heat  conduction  vector. 

In  the  work  that  follows  the  material  behavior  involves  polycrystalline  solids  that 
deform  elastically  and  viscoplastically.  Therefore,  it  is  convenient  to  decompose  the  rate 
of  deformation  tensor  additively  into  elastic  and  viscoplastic  parts 

D=De+Dp  (4.5) 

where 
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(gradv+  (gradv)T). 


(4.6) 


4.2  Elastic  Rate  of  Deformation 

The  previously  described  conservation  equations  represent  the  mathematical 
descriptions  of  stress,  strain  and  rate  of  deformation  and  they  apply  to  all  continuous 
media.  On  the  other  hand,  the  constitutive  equations  characterize  the  individual  materials 
and  their  specific  reactions  to  applied  loads.  Due  to  the  complexity  of  the  material 
response  of  different  materials  it  is  not  feasible  to  have  a  single  constitutive  relation  that 
can  be  applied  to  all  materials  under  all  loading  conditions.  For  this  reason,  the 
constitutive  model  developed  as  part  of  the  current  dissertation  work  is  tailored  to  high 
shock  compression,  high  strain,  high  strain  rate  and  high  temperature  deformation 
problems. 

In  a  general  sense,  a  constitutive  model  is  defined  as  the  function  or  functions  that 
relate  the  stress  to  the  strain  or  strain  rates.  For  the  simulation  of  hypervelocity  impacts, 
the  most  widely  used  form  of  the  constitutive  model  is  the  one  that  relates  stress  to  the 
rate  of  deformation.  In  addition,  for  the  model  presented  in  this  dissertation,  the 
constitutive  model  will  be  combined  with  thermoelastic,  thermo-visco-plastic  equations 
and  a  micro-damage  model  that  will  couple  the  elastic  moduli  and  viscosity  coefficients 
to  temperature  and  micro-damage  levels. 

The  high  compression  caused  by  hypervelocity  impact  produces  large  elastic 
volume  strain.  The  elastic  shear  deformations  however  are  small  order  because  of 
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limitations  imposed  by  onset  of  plastic  yield  in  polycrystalline  materials.  It  is 
advantageous  therefore  to  decompose  De  into  volumetric  and  deviatoric  components 

D*  =  (trDe  )l  +  D'e ,  (4.7) 

3K 

where  the  elastic  rate  of  volume  strain  is  related  to  the  mean  stress  and  the  bulk  modulus, 
K,  while  the  elastic  rate  of  small  shear  is  determined  by  the  stress  deviator  and  the  elastic 
shear  modulus,  p.  Since  the  shear  strains  must  be  small  in  the  elastic  range  it  follows 
from  generalized  Hooke’s  Law  that 

D'e  =  — T'  (4.8) 

2^ 


(v\ 

where  T'  is  the  deviator  of  the  Cauchy  stress  T,  and  ... 


Jaumann  time  rate,  whereby  for  any  tensor  A  is  given  by 


represents  the  frame  indifferent 


v  3A 

A  =  —  +  A  •  gradv  +  A  •  W  -  W  •  A  (4.9) 

dt 

where  W  is  the  rotation  tensor 

W  =^- (gradv -(grad  v)T).  (4.10) 

During  high  shock  compression  and  release,  the  large  elastic  compressive  volume 
strain  rate  is  related  to  the  applied  mean  stress  rate  and  the  elastic  bulk  modulus  K,  such 
that 

l(trD-)i=l^i 

3V  ;  3K 


(4.11) 
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where  here  and  throughout  what  follows,  o.  signtfies  the  nrean  stress  We  note 

that  in  keeping  with  the  notation  and  language  commonly  employed  in  d.scuss.ons  of 
Rankine-Hugonio,  equa.tons  introduced  in  Section  3,  the  mean  stress  there  is  referred  to 

as  -pressure'  -  P  •  ^trT .  From  expressions  (4.8)  and  (4.1 1) 

Dc  J_Hni  +  — T 

3  K  2(i 

The  pressure  associated  with  the  large  volume  stratns  tha,  accompany  high  shock 
compress, on  ,s  determined  by  the  Gruneisen  Equatton  (3.21),  here  reproduced  again  as 

Equ.  (4.13) 


(4.12) 


PM 


ll-sv]2 


l-K 
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Jo1 


V 


M' 


Yo+i 


vl-V, 


pE(v|/,T),  V>0  (4.13) 


where  il  will  be  recalled  that  1  /'po  1S  ^ 

compress, on  wave  induced  by  the  hypervelocity  impact  is  followed  by  a  reflected  tensile 
wave  from  the  targe,  rear  free  surface,  causing  dilatation  (negative  volume  strain,. 
During  the  shock  compress, on  stage  the  shock  wave  speed  is  greater  that  the  bulk  sound 
speed  (Figure  2.1,  because  of  the  large  volume  compress, on.  During  the  d, la, anon 
however,  the  crystal  lanice  structure  is  expanding,  thereby  reducing  the  tensile  wave 
speed  (Whtch  is  no  longer  a  shock  wave)  to  the  bulk  sound  speed  c„.  Thts  ,s  equivalent  to 
setting  s=0  in  Equation  (2.1).  Therefore,  during  dilatation  the  pressure  has  the  value 


P(v,T)  =  Kv 


I— yjv  (_2L|  +  a(— ! 

>  ,T  !  i-vj  li-v. 


n21 


+ 


To  +A 


1-v, 


pE(v,T).  (414) 
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The  tensile  volume  strains  cannot  assume  values  as  large  as  the  volume  strains  in 
compression  (that  can  reach  values  as  high  as  thirty  to  forty  percent),  because  of 
impending  fracture.  It  may  turn  out  that  terms  of  third  order  (and  possibly  of  second 
order)  in  vg,  in  expression  (4.14)  small  enough  to  be  negligible. 

All  polycrystalline  maierials  contain  microvoids  and  microcracks.  Since  impact 
at  hypervelocity  can  produce  large  adiabatic  temperature  increases,  the  corresponding 
dynamic  fracture  will  be  ductile,  caused  essentially  by  nucleation,  growth  and  joining 
together  of  microvoids.  The  measure  of  ductile  material  microdamage  can  therefore  be 

described  by  microvoid  volume  fraction  £  =  vy/ ,  where  vvis  the  void  volume  for  a 

material  element  having  volume  v .  The  microdamage  is  viewed  as  a  scalar-valued 
continuous  point  function  of  position  and  time,  i.e.,  £(x,t)  similar  to  the  other  field 
variables.  For  polycrystalline  metals  the  average  initial  microvoid  volume  fraction 
ranges  from  q0  ~  10~*  - 10"3  (Curran  et  al.  1987).  The  large  volume  strains  of 
compression  following  impact  ,vill  reduce  them  funner  to  insignificant  values.  On  the 
other  hand  the  dilatation  driven  by  the  tensile  mean  stress  wave,  when  high  enough, 
induces  microvoid  nucleation  and  growth,  i.e.,  material  microdamage  that  is  revealed  in 
optical  micrographs  (Curran  et  al.  1987)  and  by  the  reduced  values  for  the  elastic  moduli 
(Eftis  1996). 

In  circumstances  that  involve  hypervelocity  impact,  formulation  of  constitutive 
relations  for  the  elastic  moduli  requires  consideration  of  high  volume  strain,  large 
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temperature  changes  and  microdamage  softening.  With  this  in  mind,  the  elastic  bulk  and 
shear  moduli  can  be  represented  initially  by  the  following  form. 


K(v,E,T)=K0fK1(v)fK!(5)f„(T) 

(4.15) 

M(V,5,T)=H0f>,Wf„1fe)f„)(T) 

(4.16) 

Where  K0  =  p0c02  and  /jo  are  the  elastic  bulk  and  shear  moduli  in  the  reference  state: 
v|f=0,  T=To=300°K  and  and  where  fki,...  fk3,  and  fMi,...  fpjare  initially  undetermined 
functions  of  their  respective  arguments. 

The  bulk  modulus  is  defined  as  (Malvern  1969), 


1  do 
K  =  - — = 


p  cu 


(4.17) 


In  calculating  this  derivative  the  frequently  used  approximation  Py(h')  =  p0y0  is 
employed  (Herrmann  1969),  which  is  equivalent  to  taking  A=0  in  Equations  (3.20)  and 
(4.13)  With  this  simplification,  the  derivative  (4.17)  yields: 

(l-  y)[l+(s-y0M 


Pcct 


0-sy)3 


K0fk,(v). 


(4.18a) 


PccoO  “  v)0  -7oV)=  K0fkl(v), 


V  <  0 


(4.18b) 


Highly  different  constitutive  models  for  the  shear  modulus  for  highly  compressed 
solids  have  been  proposed  by  Steinberg  et  al.  (1980),  and  Rubin  (1987).  Alternatively, 
considering  the  lack  of  experimental  data  it  is  reasonable  to  assume  that  the  shear 
modulus  changes  with  respect  to  volume  strain  in  the  same  manner  as  the  bulk  modulus 
(Johnson  et  al.  1991  and  Lipkin  et  al.  1977),  i.e., 
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f„i(v)  =  fkiW, 


(4.19) 


where 

.  3(‘  ~  2v)  K 

2(1 +  v) 

and  vis  the  Poisson  ration. 


(4.20) 


The  effect  of  microvoid  damage  in  degrading  the  elastic  moduli  can  be  described 
by  the  well  known  hollow  sphere  proposed  Mackenzie  (1957),  and  subsequently 
modified  by  Johnson  (1981),  according  to  which 


Wi  z£j 

4/4,  +  3K0f 


U0  = 


0-0 


9Ko  +8(i, 


0  J 


(4.21) 


(4.22) 


When  '=0,  Eqs.  (4.21)  and  (4.22)  reduce  to  unity,  and  for  c=l,  (i.e.  uu  =  u)  they 
reduce  to  zero. 

The  temperature  variation  of  the  buik  and  shear  moduli  at  temperatures  above 
room  temperature  are  known  to  be  linear,  or  near  linear,  respectively,  with  the  shear 
modulus  reducing  to  zero  at  melt  temperature  (Steinberg  et  al.  1980,  Huntington  1958 
and  Lakkad  1971).  Thus 


f,K  = 


(4.23) 


^  T  T  ^ 
1  _  *o 


T  -Tr 


O  / 


(4.24) 


...I 


where  a2  and  p  are  material  constants,  T0  is  the  temperature  at  which  the  modulus  Ko  or 
p0  are  measured  and  Tm  is  the  melt  temperature  of  the  material.  The  variation  of  the  melt 
temperature  with  large  volume  reduction  can  be  calculated  according  to  the  model 
proposed  by  Steinberg  el  al.  (1980), 

Tm  =T1b0  exp(2a,v)(^)2^°"*' ,¥>0.  (4.25) 

where  Tm0  is  the  melt  temperature  at  zero  volume  strain  and  ai  is  a  material  parameter. 

The  elastic  moduli  expressions  (4.15)  and  (4.16)  can  now  be  given  in  explicit 
form,  where  for  volume  compression,  y  >  0 


„  „  ro-v)Ms-Yo)vfi  r,  „  t t  i 


r(l-vy)[l+(s-Yo)v] 
\  0-sv)3 


T  -Tn 


4HoM) 
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(4.26) 

4p0+3K0^J 

6K.+12P,  V 

(4.27) 

'  91^+8^  qJ_’ 

whereas  for  zero  volume  strain  and  dilatation,  y  <  0 ,  the  parameter  s  is  set  equal  to  zero 
in  Equations  (4.26)  and  (4.27).  For  high  shock  compression  the  microdamage 
contribution  reduces  essentially  to  unity,  while  for  small  dilatation  the  volume  strain 
makes  no  significant  contribution. 

4.3  Viscoplastic  Rate  of  Deformation 

The  plastic  rate  of  deformation  tensor  can  be  decomposed  in  terms  of  its 
volumetric  and  deviatoric  components. 


Dp  =  i(trDp)l  +  D,p 


(4.28) 
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where  -(trDp)l  represents  the  dilatational  rate  of  deformation  tensor  associated  with 
3 

void  volume  change,  1  is  the  identity  tensor,  D'p  is  the  deviatoric  rate  of  plastic 
deformation  tensor. 

Tc  define  the  volumetric  plastic  rate  of  deformation  an  increment  of  dilatational 
micro  void  volume  strain  defined  by  -  dvv  =  (l/vv  )dv  leads  to  the  logarithmic  expression. 


V,  =  "ln 


Izkl 

1-U 


(4.29) 


where  ,  Vv  represents  the  volume  strain  associated  with  void  volume  fraction  change  and 


l  =  represents  void  volume  fraction. 

It  can  be  shown  by  using  the  continuity  Equation  (4.1)  and  Equation  (4.29)  that. 


i(tr  Dp)l  =  { 
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(4.30) 


The  constitutive  model  that  describes  the  deviatoric  viscoplastic  response  of  the 
material  at  high  strain  rate  is  based  on  the  work  by  Perzyna  (1986). 


f  }' 
12--1 
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T’ , 


(4.31) 


where  the  yield  condition  requires  that 

F  =  — -1  >  0.  <«2) 

K 

The  stress  deviator  invariant  J'2  =  -^tr(T’  T'),  q  is  the  material  viscosity  and  k  is  a  yield 
function  bused  on  the  Johnson-Cook  model  (Johnson  ct  qI.  1983),  but  with  the  addition  of 
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a  damage  softening  term  to  reflect  the  deterioration  of  the  yield  stress  as  the  void  volume 
fraction  increases, 


(4.33) 


CF  is  the  equivalent  plastic  strain,  £p  the  equivalent  plastic  strain  rate,  e0  =1 .0  s'1,  Tm  is  the 
melt  temperature,  £,  the  current  void  volume  fraction,  £,o  the  initial  void  volume  fraction, 
the  critical  void  volume  fraction  at  which  the  material  locally  fractures,  and  C|,  C2.  C3, 

n  and  P  are  material  parameters.  The  term  C ,  +  C 2  (e p )"  J  describes  the  nonlinear 
isotropic  strain  hardening  characteristics  of  the  material;  term 


incorporates  strain  rate  effect;  and  term 


defines  the 


thermal  softening  characteristics  of  the  material.  When  T=Tm  the  thermal  softening  effect 


reduces  the  equivalent  yield  stress  to  zero.  The  last  term,  1  -  -f — defines  the 

damage  softening  effects  associated  with  the  presence  of  microvoids.  When  the 

damage  softening  effect  reduces  the  equivalent  yield  stress  to  zero. 

Experimental  data  is  limited  and  considerable  discrepancy  exists  in  the  reported 
values  of  material  viscosity.  Various  shock  experiments  have  shown  that  material 
viscosity  can  be  influenced  by  temperature  (Carroll  et  al.  1986),  strain  rate  (Swegle  et  al. 
1985)  and  mean  stress.  A  model  proposed  by  Carroll  et  al.  (1986)  assumed  an 
exponential  form  of  the  temperature  variation  of  viscosity  for  high  strain  rates. 
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Chhabildas  et  al.  (1979)  extracted  material  viscosity  values  for  various  materials  under 
high  shock  compression  and  elevated  temperature  from  measurements  of  shock  rise  time. 
The  model  used  in  the  present  work  assumes  that  the  lower  bound  value  for  the  viscosity 
at  melt  temperature  remains  constant  at  r)  =»  1 0*2  Poise  (Carroll  et  al.  1986), 
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(4.34b) 


where  rj0  is  the  viscosity  at  the  reference  temperature  (room  temperature)  and  high  strain 
rate  of  deformation  i.e.  1 06- 1 08  s'1,  rim  is  the  viscosity  at  melt  temperature  (Figure  4.1). 


Figure  4.1  Viscosity  as  a  function  of  temperature  for  Aluminum  applicable  to 
high  strain  rates  (ti0=1O4  P,  Tim=10'2  P  and  Tm=923  °K)  (Mineev  et  al.  1997). 
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Combining  Equations  (4.28),  (4.30)  and  (4.31),  the  plastic  rate  of  deformation 
tensor  can  be  expressed  as, 


D 


p  _ 


1 


n(T) 


J- 


-1 


T  +  [3(l-5)tl  +  ln(,-%t)j}1 


(4.35) 


This  equation  differs  from  the  one  proposed  by  Perzyna  (1986)  in  that  the 
volumetric  and  deviatoric  components  are  not  explicitly  coupled.  In  the  form  proposed 
by  Perzyna  the  constitutive  model  gives  a  non-zero  value  to  the  volumetric  contribution 
from  microvoid  growth  even  when  the  microvoid  volume  rate  of  change  is  zero.  This 
problem  is  corrected  by  using  the  form  presented  in  Equation  (4.35).  It  should  be  noted 
that  Equation  (4.35)  applies  only  for  temperatures  below  the  melt  temperature.  The 
equivalent  plastic  strain  is  defined  as  (Malvern  1969), 


(4.36) 


c'  =  J'  2lS\(n„Y  |dl,  Il0  =i[(tr  D’):  -tr(Dp  D") 

At  melt  temperature  T=Tm  and/or  at  the  critical  void  volume  fraction  it 

follows  from  (4.33)  tha*  <->0  so  that 

DPiT.i.  ->»  <4-37) 

represents  the  condition  for  local  failure,  i.e.,  the  material  loses  its  capacity  to  resist  load. 


4.4  Microdamage  Evolution 

At  the  microstructural  level,  all  polycrystalline  solids  fracture  by  the  nucleation, 
growth,  and  coalescence  of  voids  (pores)  and  cracks  (Curran  et  al.  1987).  Because  of  the 
elevated  temperatures  associated  with  hypervelocity  impact  only  the  ductile  failure  mode 
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associated  with  nucleation,  growth  and  coalescence  of  voids  (Johnson  1981)  is 
considered. 

Some  of  the  first  models  developed  to  predict  porosity  changes  of  a  material 
subjected  to  high  shock  compression  werp  developed  by  Carroll  et  al  (1972).  Their 
development  included  the  idealization  of  a  void  as  a  spherical  hole  inside  a  matrix  of 
solid  material.  From  this  initial  study,  the  model  has  been  improved  by  several  authors 
(Johnson  1981,  Carroll  1985  and  1986,  Perzyna  1986,  and  Eftis  1991,  1992  and  1996 
amongst  others).  The  development  of  the  micro-damage  model  presented  in  this 
dissertation  improves  upon  the  previous  models  in  several  ways  that  will  be  discussed  as 
the  model  is  presented. 

The  presence  of  microvoids  in  polycrystalline  materials  in  the  as-produced  state 
has  been  documented  in  many  research  and  review  articles  (Curcan  et  al.  1987).  These 
microvoids  can  vary  in  size  from  10‘-tol0'5  cm,  with  an  average  density  of  the  order  of 
106  per  cm3  and  an  average  initial  void  volume  fraction  between  1  O'3- 1 0-4  (Eftis  1996).  In 
addition,  it  has  been  demonstrated  experim  entally  that  it  is  the  growth  of  this  initial 
density  of  microvoids  in  combination  with  the  nucleation  and  growth  of  new  microvoids 
that  eventually  leads  to  the  formation  of  ductile  cracks  and  spallation  or  fracture  of  the 
material.  Figure  4.2  illustrates  the  final  effect  of  this  microdamage  process.  The 
nucleation  process  is  dependent  on  the  number  and  character  of  the  nucleation  sites  that 
are  available  in  the  crystalline  material.  In  general,  the  nucleation  of  voids  occur  at 
microscopic  heterogeneities  in  the  material  like  preexisting  flaws,  inclusions  and  second 
phase  particles,  grain  boundaries  and  subgrain  structures  (Curran  et  al.  1987).  Figure  4.3 
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is  evidence  that  inclusions  in  the  material  can  serve  as  nucleation  sites.  In  this  figure  the 
inclusions  are  still  observable  inside  the  voids  that  have  grown  to  many  times  their  initial 
diameters. 


Figure  4.2  Ductile  cracks,  (a)  Ductile  crack  propagation  by  void  coalescence;  (b) 
Tip  of  ductile  crack  shown  in  (a)  at  higher  magnification  (Curran  et  al.  1987). 


Figure  4.3  Fracture  surface  of  ETP  Copper  (Zurek  et  al.  1996). 


4.4.1  Microvoid  Nucleation  Rate 

The  model  used  for  the  simulation  of  void  nucleation  has  its  origin  in  the  form 
proposed  by  Seaman  et  al.  (1976).  Perzyna  (1986)  further  modified  this  model  to  an 
activation-rate  form,  which  also  makes  provision  for  the  effect  of  void  interaction  during 
the  nucleation  process,  so  that  the  void  nucleation  rate  model  has  the  form 


In  the  above,  m:io-on|  is  the  nucleation  activation  energy,  am  is  the  mean  stress.  on  is 

the  threshold  mean  stress  for  nucleation,  T  is  the  temperature,  k  is  the  Boltzman  constant 
and  m;  are  material  parameters  and  h  simulates  the  effect  of  void  interaction.  The  void 
nucleation  contribution  to  the  total  void  volume  fraction,  given  by  the  time  integration  of 
Equation  (4.38),  has  a  limiting  value,  i.e.,  t,n  <£,nm.  This  limit  is  controlled  by  the 


number  of  nucleation  sites  available  in  the  material  in  the  form  of  heterogeneities  (Curran 
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el  al.  1987).  The  values  for  the  nucleation  threshold  stress  found  in  high  velocity  impact 
literature  are  mostly  obtained  from  high  velocity  plate  impact  tests  (Curran  et  al.  1987). 
The  high  shock  compression  associated  with  impact  at  hypervelocities  can  cause 
substantial  strain  hardening  of  the  material,  particularly  at  the  impact  sides  of  the  target 
and  projectile.  This  effect  is  indirectly  exhibited  by  the  higher  nucleation  thresholds  for 
the  harder  metals,  e.g.,  Ni-Cr-Steel,  as  compared  to  softer  materials  such  as  Copper  or 
Aluminum,  and  may  be  expressed  by  assuming  the  microvoid  threshold  nucleation  stress 
to  be  proportional  to  the  current  yield  stress, 


where  cn0  is  the  nucleation  threshold  determined  by  plate  impact  tests  and  the  rest  of  the 
parameters  are  as  defined  in  Equation  (4.33).  The  nucleation  threshold  has  the  same 
functional  form  as  the  Johnson-Cook  model  for  the  equivalent  yield  stress,  except  that  the 
temperature  factor  is  excluded  since  the  temperature  dependence  is  already  included  in 
the  void  nucleation  rate  equation  ?vove.  Several  authors  have  implemented  this  model 
for  nucleation  but  have  used  a  constant  nucleation  threshold  (Perzvna  1986,  Eftis  et  al. 
1991  and  1992,  and  Eftis  1996).  Since  most  of  the  validation  of  the  model  was  done  with 
simulations  of  high  velocity  plate  impact  tests,  the  results  that  were  obtained  were 
acceptable.  However,  when  a  constant  value  for  the  threshold  is  used  to  simulate 
hypervelocity  impact  the  nucleation  model  with  constant  nucleation  threshold  stress  no 
longer  gives  acceptable  results,  and  a  variable  threshold  as  described  by  Equation  (4.39), 
is  necessary. 
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4.4.2  Micro  Void  Growth  Rate 

The  rate  of  change  of  the  void  volume  fraction  model  is  based  on  the  idealized 
spherical  hole  described  by  Carroll  et  al.  (1972).  In  their  model  they  assumed  that  the 
yield  stress  of  the  material  was  constant  (perfectly  plastic  material)  and  ignored  thermal 
effects.  The  spherical  void  model  was  further  developed  by  Johnson  (1981)  to  include 
strain  rate  effects  in  the  constitutive  model  of  the  material  and  was  the  first  to  infer  that 
thermal  effects  might  be  important.  The  model  was  later  improved  by  Carroll  by 
studying  the  effects  of  several  material  isotropic  hardening  laws.  Perzyna  (1986)  included 
linear  strain  hardening  in  addition  to  the  rate  effect.  The  linear  hardening  model 
proposed  by  Perzyna  was  proved  unsatisfactory  by  Eftis  et  al.  (1991,  1992)  and  was 
replaced  by  a  saturation  hardening  model.  Thermal  effects  were  not  included. 

The  development  of  the  spherical  void  model  is  presented  and  includes  several 
improvements  over  the  models  presented  before.  The  derivation  of  the  model  begins  by 
replacing  the  microvoids  embedded  in  the  material  by  a  random  distribution  of  spherical 
voids  (Figure  4.4). 


-a 


Figure  4.4  Element  of  the  solid  with  microvoids  replaced  by  idealized  spherical 
voids  (Eftis  1996). 
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If  the  rectangular  volume  of  Figure  4.4  is  subjected  to  the  actions  of  mean  stress 
a(t)  that  varies  with  time,  only  the  solid  cross-sectional  area  will  support  the  applied 
stress.  The  sectional  area  occupied  by  the  voids  does  not  support  any  stress.  Thus. 
mecKnical  equilibrium  requires  that  (from  this  point  on  the  time  dependence  of  the 
variables  is  implicit), 

Ao  =  Aa  (4.40) 

where  A  is  the  total  cross-sectional  area  (i.e.,  including  the  area  of  the  microvoids),  a  is 
the  mean  stress  acting  on  A ,  A  is  the  solid  portion  of  the  cross-sectional  and  o  is  the 
mean  stress  acting  on  the  solid  material  with  sectional  area  A.  For  a  random  distribution 
of  voids  shapes  and  sizes  the  following  holds, 

A  v 

A=?  (4-41> 

where  v  is  the  volume  of  solid  material  and  v  is  the  total  volume  of  material. 

Substituting  Equation  (4.41)  in  (4.40), 

a  =  aa  (4.42) 

where,  a  =  V  is  defined  as  the  void  distension. 

Considering  an  idealized  void  with  spherical  shape  and  radius  a  in  a  sphere  of 
radius  b  where  the  outer  boundary  is  subjected  to  mean  stress  cr  (figure  4.5). 
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Figure  4.5  Idealized  spherical  void  with  average  mean  stress. 


The  equation  of  motion  for  the  configuration  shown  in  Figure  4.5  can  be  defined 
in  spherical  coordinates.  Taking  account  of  spherical  symmetry, 

<3T  2  (T  r  \  •• 

—  +  “(T„  -Tee)  =  pr  (4.43a) 

ct  r 

The  boundary  conditions  require 


T„(m)-0 

T„(b,t)-_o 


^  (4.43b) 


where  T„  is  the  radial  stress  acting  on  a  spherical  differential  element  (Figure  4.5).  Too  is 
the  tangential  stress,  p  is  the  density  of  the  solid  material,  and  r  is  the  material  radial 
acceleration.  Making  the  assumption  that  nearly  all  the  volume  change  is  associated  with 
void  growth  and  very'  little  with  volume  strain  in  the  solid  material  surrounding  the  void, 

r’  =ro3  -B(t)  (4.44) 

b3 

b 3  —  a 3 


G  = 


(4.45) 


52 


where  B(t)  is  a  function  related  to  the  rate  of  void  growth. 

By  means  of  Equation  (4.44)  the  radius  of  any  point  that  initially  had  radius  ro, 
can  be  determined.  Differentiating  Equation  (4.44)  twice  with  respect  to  time. 


f.J&.ME.  -  ( 

3r2  9r5 

Integrating  Equation  (4.43)  from  a  to  b  (inner  to  outer  current  radius)  and  substituting 
Equation  (4.46), 


p[4J(b,  t)-vf'(a,  t)]=  ao  + 


(4.47a) 


where. 


3r  18r 


(4.47b) 


By  means  of  Equations  (4.44)  and  (4.45), 


B(t)  =  a0J 


Qq  -a 

Lao  “I  . 


(4.48) 


b  3 

where  ao  is  the  initial  inner  radius  and  a0  =  — p5 — -  the  initial  void  distension. 

bo  -ao 

Substituting  Equation  (4.48)  into  (4.47)  and  rearranging  terms, 

-  br2, 


TQ(d.a,a)  =  ao+  f-(T„  -Tee)dr 
:  r 


(4.49a) 
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(4.49b) 


Q(d,a,a)=  [(a  —  l)-10  -a~I/3]a-  — [(a-l)~4/3  -a~4/3]d:.  (4.49c) 

6 

At  this  point,  a  constitutive  model  for  the  material  surrounding  the  void  has  to  be 
introduced  to  completely  define  Equation  (4.49).  It  can  be  shown  that  the  elastic  and 
elastic-plastic  phases  of  the  deformation  of  the  material  surrounding  the  void  produce  a 
very  small  volume  change  compared  to  the  plastic  deformation  contribution.  For  this 
reason,  as  an  approximation,  void  growth  is  going  to  be  modeled  as  a  fully  plastic  process 
(Johnson  1981).  In  other  words,  void  growth  is  only  going  to  occur  when  all  the  gray 
area  shown  on  Figure  4.5  has  yielded.  The  yield  function  Y  used  for  the  constitutive 
model  has  a  form  similar  to  the  modified  Johnson-Cook  model  shown  in  Equation  (4.33) 
but  w  ithout  the  damage  softening  effect  and  w'ith  the  '•train  rate  effect  as  an  additive  term. 

Y.(c'.c'.T.r)=--!r(T„-Ts)  =  ±[ci+C!(E')n]l-fihji)  ±n(T(r))e'  (4,50) 

V-5  ‘m  *0  9 

where  rj(T ( r) )  is  the  temperature  dependant  viscosity  at  any  radius  of  the  material 
surrounding  the  spherical  void.  The  negative  sign  in  the  difference  of  stresses  term 
comes  from  the  fact  that  when  the  applied  mean  stress  (a)  is  positive,  the  difference 
between  the  radial  and  tangential  stresses  is  negative  (Malvern  1969).  The  ±  sign 
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indicates  that  the  yield  stress  is  positive  when  the  mean  stress  is  positive  and  negative 
when  the  mean  stress  is  negative. 

In  applications  of  the  microvoid  growth  model  for  description  of  spall  fracture 
induced  by  plate  impact  (Johnson  1981  and  Eftis  et  al.  1991)  it  ;  found  that  the 
numerical  simulations  w'ould  correlate  with  experimental  data  only  it  the  viscosity  of  the 
solid  material  surrounding  the  void  has  a  value  of  approximately  10  Poise,  three  to  four 
orders  of  magnitude  smaller  than  room  temperature  values.  It  was  suggested  that  the 
large  shear  strains  at  the  pore  wall  cause  temperature  increase  close  to  melt  temperatures, 
thereby  greatly  reducing  the  viscosity  of  the  surrounding  material.  The  variation  of  the 
viscosity  with  temperature  is  given  by  Equation  (4.34),  with  the  melt  temperature 
Tm=Tm0,  i.e.,  the  value  at  zero  volume  strain  (Equation  4.44). 

To  simplify  the  integration  of  Equation  (4.50)  when  substituted  into  (4.49a).  the  absolute 

temperature  used  in  Equation  (4.50)  is  taken  at  a  fixed  radius  r(t)  in  the  vicinity  of  the 
inner  radius  (a), 

r(t)  =  a -q(b-a).  0<q<l.  (451, 

If  q=0  then  r(t)  is  the  inner  radius  (a),  and  if  q-1  then  r(t)  represents  the  outer  radius 
(b). 

Substitution  of  Equation  (4.51)  into  (4.50)  eliminates  the  radial  dependence  of  the 
temperature  and  implicitly  of  the  viscosity  as  well.  Thus 


jr(T„  ~tm)=±[c,  +C.(cp)" 


1- 


f  T-T  ^ 
1  10 

“To  J 


±q(T)ep  (4.52) 


where  T  is  the  temperature  at  radius  r. 
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From  Equation  (4.44)  the  radial  displacement  can  be  calculated  to  be 


and  for  small  displacements  and  strains,  the  small  strain  tensor  in  spherical  coordinates 
and  with  radial  symmetry  has  the  matrix  representation 
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The  equivalent  plastic  strain  is  given  by  (Malvern  1969), 

c'=^|(E;-Ei)1+(E;-E;)ir 

J  ' 

which  reduces  to 


(4.54) 


(4.55a) 


3  :  r  cr 


The  equivalent  plastic  strain  rate  is 


(4.55b) 


er  = 


f2k. 

au,  'i 

Ui  r 

CT  J 

(4.56) 


The  integral  that  appears  on  the  right  hand  side  of  Equation  (4.49a).  after 
substitution  of  Equations  (4.55),  (4.56)  and  (4.53)  in  Equation  (4.52),  is  not  developed 
here  because  it  is  a  lengthy  calculation.  Only  the  Final  result  is  shown. 
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“  F(a1  (4.57a) 


where 


B(t)/b’  . 

I(a)=  J  -|xO  +  x)"M|"dx 

B(t )  a’  ^ 


(4.57b) 


f (t)  =  i-f  T  T°- 


T  -T0 


(4.57c) 


F  =  a 


a0  -1 
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(4.57d) 


The  integral  in  Equation  (4.57b)  is  non-elementary  and  has  to  be  solved  numerically. 

The  substitution  of  Equation  (4.57)  in  Equation  (4.49)  gives, 

iQfed.a)=(a6)+^  c,lr(^j-(^  C;I(a)  f(T)+^n(T)^jF(a)  (4.58) 


If  we  take  the  limit  of  Equation  (4.58)  as  a  and  a  approach  zero. 


IcJ  3 


2V3  rif  a 
-  C.ln  - 


I  C,I(o)  f(t) 


(4.59) 


Which  represents  the  threshold  mean  stress  for  void  growth  (+),  or  void  volume  reduction 
(-),  and  includes  the  power  isotropic  hardening  and  the  thermal  softening.  Once  the 
externally  applied  mean  stress  (ct) reaches  the  threshold  (ctg)  the  spherical  void  begins 
to  grow  or  shrink. 

Substituting  Equation  (4.59)  into  (4.57), 


xQ(d,d,a)  =  a(o  -  oc)T  — j=t\(t 


4  f—i  a  -/ 


3-v/3  V  ]  a(a  -  l) 


(4.60) 
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To  express  Equation  (4.60)  in  terms  of  the  void  volume  fraction  (£,)  instead  of 
void  distension  (a),  the  following  relation  can  be  derived  from  the  definitions  of  t  and  a 
mentioned  above, 


1 

a  = - 

1  t 


Differentiating, 


2*1  i 

mT +o-?)! 


Then  substituting  Equations  (4.61)  into  (4.60)  and  arranging  terms, 


iM  +  B^itlMlFOOUto-OG) 


where. 


(4.61a) 

(4.61b) 

(4.61c) 

(4.62) 

(4.63) 

(4.64) 

(4.65) 

(4.66) 
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(4.68a) 


(4.68b) 


(4.68c) 
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The  term  t(a,^  +  B,42)  represents  the  inertial  resistance  to  microvoid  volume 
change,  r|(T ):^F(^)|  the  viscous  material  resistance  to  plastic  flow  and  (a  -  oG )  the  mean 


stress  difference  that  drives  the  void  growth. 

As  voids  grow  in  the  material,  they  initially  grow  independently.  However,  in  the 
later  stages  before  fracture,  they  begin  to  interact.  First,  their  plastic  strain  fields  intersect 
and  then  a  stronger  interaction  takes  place  (Curran  et  al.  1987).  It  is  this  interaction  that 
represents  the  final  stage  of  failure  of  the  material.  Eftis  et  al.  (1991)  proposed  a 
modification  to  the  void  growlh  Equation  that  would  take  void  interaction  into  account  as 
the  growth  process  progresses.  This  is  done  by  the  inclusion  of  an  exponential  function 
in  Equation  (4.62), 

*(A,q  +  B,42)±q(T)|4F(^)|  =  e*(o-60),  ^>0  (4.69a) 

t(a,  \  +  B,42  )±  rj(T =  (o  -  oG ) ,  S<0 


(4.69b) 
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where  (p  is  a  material  parameter.  Equation  (4.69b)  differs  from  (4.69a)  in  that  it  does  not 
include  the  coalescence  effect  since  under  void  volume  reduction  there  is  no  coalescence 
effect  of  the  microvoids. 

The  microvoid  nucltation  rate  given  by  Equation.  (4.38),  when  added  to  the 
solution  of  the  Equation  (4.69)  for  the  void  growth  rate  4g,  provides  the  total  microvoid 


4.4.3  Temperature  Near  the  Wall  of  a  Microvoid 

The  temperature  change  in  the  vicinity  of  the  walls  of  a  microvoid  (hollow 
sphere)  is  governed  by  the  plastic  rate  of  mechanical  work, 

pCJ  =  cotr(T-Dp/  (4.71) 

where  C,  is  the  specific  heat  at  constant  volume,  and  w  gives  the  percentage  of  the  rate  of 
plastic  work  that  is  not  stores  as  cold  work,  i.e.,  that  dissipates  into  heat.  For 
polycrystalline  materials  coM).9.  In  view  of  the  spherical  symmetry  only  the  radial 
component  of  T'  and  Dp  are  non-zero  so  that 

tr(T'Dp)=T;-Dp„  (4.72) 

where 
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d--!% 

3  r3 


(4.73) 

(4.74) 


and  by  using  the  simplification  proposed  by  Equation  (4.51)  it  can  be  shown  that ' 


(Tn-Tj  =  ±V3[ci+C2(ep)n]  1 


T-T0 

Tm-T0 


±ti(t)^p 


(4.75) 


Further  omitting  the  detailed  calculations  it  can  be  shown  that 


(4.76) 


t'-iM  ir1+®ar  +ifi+5fc)' 

3  r3  3  r3  3  r3 


(4.77) 


B(t)_ 

o 

1 

-  3 

r 

J ^ 

t 

o 

-G^ 

B(t) 

c  1 

f3 

$0-  0<p(5) 

(4.78) 


(4.79) 


tp(c)  =  <  1  +  q  -1 


(4.80) 


From  Equations  (4.72)-(4.74)  expression  (4.71)  can  be  reduced  to 


pC.T  =  coi(Tn-TM)^l 

9  r 


(4.81) 
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Equation  (4.81)  represents  the  temperature  rate  at  any  radius  that  is  a  fraction  (q). 
between  the  inner  and  outer  radius  of  the  void.  The  temperature  calculated  through  the 
integration  of  Equation  (4.81)  gives  the  temperature  used  in  Equation  (4.69a)  for  the 
temperature  dependant  viscosity  and  in  the  temperature-softenin':,  factor  (Equation  4.57c) 
applied  to  the  growth  threshold  mean-stress  in  Equation  (4.67).  This  void  temperature 
calculation  provides  analytical  justification  for  the  reduced  temperature-independent 
viscosity  used  by  Johnson  (1981)  and  Eftis  et  al.  (1991),  and  represents  an  important 
improvement  over  the  previous  models. 

Combining  expressions  (4.75)  to  (4.81)  the  temperature  rate  at  some  location  in 
the  vicinity  of  the  inner  pore  wall  can  be  expressed  by 


T  =  -^-a/J T: 

9pCvV 


*0-e> 


1  +  q 


fl  \'n 


where 


(4.82) 


n  =i(T„-T*)!  (4.83) 

represents  the  second  invariant  of  the  stress  deviator  tensor,  the  equivalent  of  which  can 
be  calculated  by  using  Equations  (4.75)-(4.79),  and  where  0<q<l.  Equation  (4.82) 
makes  explicit  the  fact  that  the  thermal  heating  near  the  wall  of  the  microvoid  is  caused 
by  the  shear  stress  and  plastic  rate  of  change  of  the  microvoid  volume  fraction. 


4.5  Equations  for  Internal  Energy  and  Temperature  Rates 

Due  to  the  rapidity  of  the  deformation  following  hypervelocity  impact,  the 
deformation  and  fracture  processes  develop  under  adiabatic  conditions.  In  addition,  if  we 
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consider  that  there  are  no  external  heat  sources,  then  the  local  energy  rate  balance 


Equation  (4.4)  simplifies  to 
pE  =  tr(TD) 


(4.84) 


This  exr  session  can  be  developed  further  by  decomposing  the  stress  and  rate  of 
deformation  tensors  into  deviatoric  and  mean  (or  volumetric)  components,  making  use  of 
(4.5)  and  the  relation 


trD'  =  — 


(4.85) 


It  the  follows  that 


Pe  =  -o„-^+tr(T-D")+tr(T.D'). 

l-v 


(4.86) 


where  the  mean  stress  is  given  by  Equations  (4.13)  and  (4.14).  If  we  specifically 
consider  the  internal  energy  as  thermal  energy,  then  pE  =  pC,T  .  Since  small  elastic 

shear  deformation  does  not  cause  temperature  change,  tr(T-D")=0.  In  as  much  as 
elastic  volume  reduction  increases  temperature  while  elastic  dilatation  causes  temperature 
reduction,  consistency  requires  -am  to  be  replaced  by  |am|.  Taking  account  of  the  fact 
that  onlv  a  fraction,  co,  of  the  plastic  work  dissipates  into  heat.  Where  o=0.9  for 


polycrystalline  metals,  we  have  the  local  temperature  rate. 


pC.T  =  |a„|-^-  +  o>tr(T-Dp), 
l-V 


(4.87) 


CHAPTER  5 


NUMERICAL  IMPLEMENTATION  OF  THE  CONSTITUTIVE  MODEL 

The  set  of  coupled,  highly  non  linear  equations  described  in  chapters  two.  three 
and  four  cannci  be  solved  analytically,  unless  simplifications  of  the  constitutive  model 
are  introduced.  For  this  reason  numerical  solution  of  the  conservation  equations  and 
constitutive  model  represents  the  only  viable  approach. 

Historically,  the  computer  programs  used  in  the  simulation  of  shock  deformation 
problems  have  been  called  “Hydrocodes”.  This  name  was  used  because  in  the  initial 
development  these  computer  codes  used  a  hydrodynamic  material  model  to  simulate  solid 
materials,  i.e.,  where  the  stress  and  strain  fields  are  assumed  to  have  no  deviatoric 
components  (Johnson  et  al.  1987).  One  of  the  justifications  for  this  type  of  simplification 
was  the  fact  that  the  shock  stresses  are  considerably  larger  than  the  dynamic  yield  stress 
of  the  material.  After  recognition  that  deviatoric  effects  were  important  in  the  simulation 
of  initiation  and  e'-pansion  of  the  detonation  products  of  high  explosives,  Wilkins  et  al. 
(1974)  developed  the  procedure  to  incorporate  deviatoric  stress  effects  into  computer 
codes.  Since  then,  a  great  volume  of  work  has  been  done  directed  towards  the 
improvement  of  the  numerical  schemes  used  to  solve  the  differential  equations  that  define 
shock-induced  deformation  problems. 

Currently  there  are  several  “Hydrocodes”  produced  by  government  and  private 

organizations.  The  codes  most  widely  used  for  the  simulation  of  hypervelocity  impact 

problems  are,  Dyna2D  and  3D  produced  by  the  Lawrence  Livermore  National 
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Laboratory’,  CTH  produced  by  Sandia  National  Laboratory',  EPIC-2  and  3  and 
Autodyn2D™  and  3D™  produced  by  Century  Dynamics  Inc. 

For  the  purpose  of  selecting  a  “Hydrocode”  to  implement  the  constitutive  model, 
evaluations  of  CTH,  Dyna2D  and  Autodvn2D™  were  performed  b~',ed  on  their 
capabilities  and  ease  of  implementation  of  new  constitutive  models.  Of  the  three  codes. 
Autodyn2D™  was  selected  because  it  has  a  user-friendly  interface,  includes  the 
capability  of  running  simulations  with  Lagrangian  and/or  Eulerian  grids,  and  has  the  built 
in  capability  of  compiling  external  functions  to  implement  new  constitutive  models. 
Even  though  the  capability  to  include  external  functions  was  an  important  factor  in  the 
selection  of  the  code,  in  the  end,  that  capability  was  not  used  because  the  constitutive 
model  described  in  this  dissertation  is  too  complicated  to  be  implemented  just  through 
external  functions.  Century'  Dynamics  Inc.,  the  originators  of  Autodyn,  provided  access 
to  some  of  the  internal  subroutines  that  needed  to  be  modified  to  successfully  implement 
the  model.  Because  of  copyright  restrictions,  a  thorough  description  of  the  subroutines 
that  were  modified  cannot  be  published.  - 

Although  the  equations  described  in  the  previous  chapters  describe  a  continuum, 
the  computer  code  has  to  use  a  discretized  version  of  the  equations.  This  discretization  is 
applied  both  in  space  and  in  time.  The  finite  element  method  (FEM)  and  the  finite 
difference  method  (FDM)  are  the  methods  of  discretization  in  space.  In  “Hydrocodes” 
the  time  integration  is  done  through  an  explicit  central  difference  scheme  in  which  the 
time  derivatives  are  replaced  by  difference  equations  where  the  value  of  the  function  at 
the  new  time  step  is  only  a  function  of  the  value  in  the  previous  time  step. 
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Once  the  discretized  form  of  the  equations  is  obtained,  they  have  to  be  structured 
within  the  framework  of  a  computer  code.  Most  of  the  codes  available  today  are  based 
on  the  work  by  Wilkins  et  al.  (1974),  with  minor  modifications.  He  described  how  the 
discretized  form  of  the  equation'  had  to  be  obtained  and  the  numerical  scheme  that  had  to 
be  used  to  advance  the  simultaneous  solution  of  the  equations  in  time  at  every  time 
increment. 

5.1  Numerical  Simplifications 

Because  the  solution  of  the  set  of  field  and  constitutive  equations  is  framed 
around  a  pre-established  numerical  scheme  that  is  the  foundation  of  every  computer  code, 
the  inclusion  of  a  given  set  of  constitutive  equation  has  to  adapt  to  the  established 
scheme,  ll  would  be  desirable  to  be  able  to  implement  any  constitutive  model  with  no 
restrictions,  but  that  is  not  possible. 

During  the  implementation  and  testing  of  the  constitutive  equations  in  the 
Autodyn™  code,  a  series  of  difficulties  were  encountered.  These  difficulties  restricted 
the  implementation  of  the  equations  and  a  series  of  modifications  had  to  be  applied  to 
make  the  simulations  feasible. 

5.1.1  Modifications  to  the  Viscoplastic  Constitutive  Model 

The  viscoplastic  constitutive  model  presented  in  Section  4.3  of  the  dissertation 
was  not  implemented  in  the  original  form  because  it  required  modifications  of  sections  of 
the  code  to  which  we  did  not  have  access.  For  this  reason,  the  calculation  of  the 
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equivalent  plastic  strain  as  described  in  Equation  (4.36)  was  not  used  and  the  calculation 
already  implemented  in  Autodyn™  was  used  with  some  modifications  (Autodyn™. 
Theory  Manual). 

T'  c  viscoplastic  constitutive  model  in  Autodyn™  takes  the  form  described  by 
Wilkins  et  al.  (1974)  with  the  additional  consideration  that  the  model  has  to  take  into 
account  the  volume  compressibility  of  the  material  due  to  the  presence  of  microvoids  that 
can  grow  or  reduce  in  size. 

The  constitutive  model  that  describes  the  deviatoric  viscoplastic  response  of  the 
material  is  based  on  the  von  Mises  yield  criteria,  which  can  be  stated  as  follows. 


a 


where  J;  is  the  second  invariant  of  the  stress  deviator  tensor  and  Y  is  the  equivalent  yield 
stress.  Thus,  when  J:  <  Y /  the  material  is  in  an  elastic  state  of  stress,  and  when 

J :  >  Y ' /.  the  material  is  in  a  state  of  plastic  yielding. 

The  definition  of  the  equivalent  yield  stress  (Y)  is  based  on  the  Johnson-Cook 
model  (Johnson  cl  a!  1983),  but  with  the  addition  of  a  damage  softening  term  to  refect 
the  deterioration  of  the  yield  stress  as  the  void  volume  fraction  increases,  i.e.. 
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which  is  identical  to  Equation  (4.33). 

The  equivalent  plastic  shear  strain  is  defined  as. 
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3p 

where,  p.  is  the  shear  modulus  as  defined  in  Section  4.2. 

It  should  be  noted  that  the  definition  of  the  equivalent  plastic  strain  defined  by  Equation 
(5.3)  has  all  the  material  behavior  features  (strain  hardening,  rate  effects,  thermal 
softening  and  material  damage  softening)  that  Equation  (4.36)  presents. 


5.1.2  Modifications  to  the  Void  Growth  Model 

After  implementation  of  the  void  growth  model  described  by  Equations  (4.62)- 
(4.68).  a  series  of  tests  were  performed  to  determine  what  would  be  the  time  increment 
that  would  give  a  converging  solution  given  the  explicit  integration  scheme  used.  It  was 
discovered  that  the  time  step  required  w'as  tw'o  orders  of  magnitude  smaller  than  the 
required  to  do  the  simulation  without  the  growth  model.  This  implied  that  a  simulation  of 
a  hypervelocity  impact  event  as  described  in  Chapter  six  would  require  approximately  30 
days.  This  was  prohibitive  given  the  time  constraints.  For  this  reason,  the  model 
presented  in  Equations  (4.62)-(4.68)  was  simplified.  The  simplification  consists  of  the 
elimination  of  the  inertial  effects,  which  reduces  the  nonlinear  second  order  differential 
equation  to  a  first  order  nonlinear  differential  equation. 


(5.4a) 
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The  implications  of  this  simplification  in  a  hypervelocity  impact  problem  are  hard 
to  infer  and  will  be  material  of  further  research.  Nonetheless,  the  modifications  and 
improvements  made  to  the  previous  models  presented  by  Johnson  and  Eftis  are  still  part 
of  the  simplified  model.  The  temperature  local  to  the  microvoid  wall,  the  temperature 
dependence  of  the  viscosity  and  the  use  of  a  viscoplastic  constitutive  equation  that 
includes  nonlinear  strain  hardening,  strain-rate  hardening  and  temperature  softening,  for 
the  behavior  of  the  solid  material  of  the  hollow  sphere  are  included  in  the  simplified  void 
growth  rate  equation  given  by  (5.4). 

The  question  of  the  importance  of  maintaining  the  inertial  term  of  the  differential 
Equation  (4.62)  is  simulations  of  very  high  strain-rate  dynamic  microvoid  growth  is  an 
open  question  at  this  time,  w-ith  adherents  on  both  sides  of  the  issue  (Johnson  1981  and 
Tong  et  a!  1995). 

5.1.3  Modifications  to  the  Temperature  Rate  Equation 

Even  though  the  temperature  raf'  Equation  (4.87)  was  derived  through  the  use  of 
the  energy  rate  balance  equation  that  applies  to  any  material  under  any  loading  condition, 
its  implementation  presented  some  difficulties.  Under  large  volume  expansion.  Equation 
(4.84)  has  the  tendency  to  give  temperatures  that  are  well  below  the  room  temperature 
and  close  to  absolute  zero.  This  was  observed  specifically  in  the  ejecta  rim  produced  in 
the  simulation  of  hypervelocity  impacts  of  spherical  projectiles  described  below.  It  is 
suspected  that  Equation  (4.87)  may  need  a  time  step  that  is  much  smaller  then  the  one 
used  for  it  to  converge  to  a  solution.  A  study  of  this  problem  will  be  material  for  further 
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research.  Thus,  the  temperature  calculation  was  done  using  the  model  already 

implemented  in  Autodyn™  (Autodyn™,  Theory'  Manual),  which  is  derived  usins 
classical  thermodynamics. 


T  = 


^LT£iT++co-!ifeA^ 

1-V  P  pCv 


(5.5) 


where  y  is  the  Gruneisen  coefficient  as  defined  by  Equation  (3.20). 


CHAPTER  6 


NUMERICAL  SIMULATIONS 


In  the  literature  review  presented  at  the  beginning  of  this  dissertation,  several 
papers  were  cited  that  employed  similar  constitutive  models  to  the  one  presented  in  this 
work  for  the  simulation  of  plate  impact  problems.  Even  though  those  models  were 
simpler  than  the  one  described  here,  the  use  of  those  models  for  the  simulation  of  impact 
problems  were  limited  to  very  simple  configurations.  Most  of  the  simulations  were  done 
for  one-dimensional  shock  plate  impact  problems  in  the  low  to  medium  velocity  range 
(Eftis  et  al  1991,  Cortes  et  al.  1995).  Simulations  of  impact  in  the  hypervelocity  range 
have  used  a  very  simple  limit  strain  or  stress  fracture  model  that  does  not  model  the 
evolution  of  material  degradation  in  a  realistic  form  (Quinones  et  al.  1998).  For  this 
reason,  the  numerical  simulation  of  hypervelocity  impact  of  spherical  projectiles  with 
metallic  target  plates  using  the  microdamage  constitutive  model  shown  here  constitutes  a 
step  forward  in  the  constitutive  modeling  and  simulation  of  hyper/elocity  impact. 

6.1  A  Shock  Pressure  and  Temperature  Calculation 

Prior  to  the  implementation  of  the  constitutive  microdamage  model  into 
Autodvn™,  a  series  of  simulations  of  a  representative  point  behind  a  shock  wave  were 
performed.  A  special  program  was  developed  using  Delphi™  (Pascal  compiler)  to  time 
integrate  the  governing  equations,  and  plots  of  different  variables  were  generated  through 


a  Matlab™  program.  The  objective  of  these  simulations  was  to  verify  that  all  the 
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components  of  the  constitutive  model  were  showing  the  expected  trends  in  terms  of  the 
behavior  of  the  field  variables.  A  comparison  is  shown  in  Figure  6.1  of  the  calculated 
and  experimental  shock-pressure  volume  strain  curves  for  Copper.  This  figure  shows  a 
very  good  agreement  between  the  calculated  and  experimental  curves  indicating  the 
suitability  of  the  equation  of  state  presented  in  Chapter  3  for  pressures  up  to  142  GPa. 
Figure  6.2  illustrates  the  calculated  temperature  rise  with  increasing  pressure,  where  the 
solid  curve  is  determined  from  the  temperature  rate  Equation  (4.87),  while  the  open  circle 
curve  was  calculated  by  others  using  the  equations  of  classical  thermodynamics 
(McQueen  ct  al.  1960).  In  addition.  Figure  6.2  includes  a  plot  of  the  corresponding  melt 
temperature  (Equation  4.25,  dashed  line)  for  the  corresponding  pressure  rise,  which 
indicates  that  at  approximately  155  GPa,  the  material  reaches  the  melt  point  according  to 
equation  (4.87)  (solid  line). 


Figure  6.1  Calculated  and  experimental  (open  circles)  pressure  versus  volume 
strain  curves  for  Copper  (circles,  McQueen  et  al.  1960). 


Figure  6.2  Calculations  of  temperature  versus  shock  pressure  for  Copper  using 
thermal  energy  rate  balance  (solid  line)  and  classical  thermodynamics  (circles,  McQueen 
et  al  1960),  and  melt  temperature  (dashed  lines). 

Because  of  the  great  difficulty  in  obtaining  temperature  measurements,  directly  or 
indirectly,  of  thermal  changes  that  are  of  nano-second  duration,  there  is  no  experimental 
udta  available  that  measures  high  shock-pressure  induced  temperature  fluctuations.  For 
this  reason,  a  conclusion  as  to  which  temperature  calculation  gives  a  temperature  closer 
to  reality  cannot  be  reached.  This  topic  requires  further  research. 


6.2  Hypervelocity  Impact  of  Soda-Lime  Glass  Projectiles  on  Alnoo  Target  Plates 
To  test  the  capability  of  the  constitutive  model  to  simulate  the  deformation 
process,  micro-damage  evolution  and  finally  fracture  developed  under  hypervelocity 
impact,  for  different  projectile  diameter  to  plate  thickness  ratios  a  series  of  experimental 


results  on  1100  Aluminum  plates  were  selected.  Horz  et  al.  (1995)  presented  results  from 
a  series  of  hypervelocity  impact  experiments  of  soda-lime  glass  spherical  projectiles  of 
3.175  mm  diameter  that  impact  Alnoo  plates  of  varied  thickness  at  a  nominal  velocity  of 
6.0  km's.  In  all  ten  experiments  were  presented  where  the  target  plate  thickness  varied 
from  1.06  to  51.75  mm,  while  all  the  other  parameters  (projectile  and  target  materials, 
projectile  diameter  and  velocity)  were  kept  fixed.  Out  of  these  ten  experiments,  three 
were  selected  to  be  simulated  using  the  Autodyn™  code  that  was  modified  to  include  the 
constitutive  model  presented  in  this  dissertation.  Figure  6.3  shows  postmortem  pictures 
of  the  three  cases  selected  from  Horz  et  al.  1995.  Also,  Table  6.1  lists  the  dimensional 
characteristics  of  the  three  hypervelocity-impact  cases.  These  three  hypervelocity 
experiments  were  selected  because  they  cover  a  wide  range  of  plate  thickness  that  result 
in  three  widely  varied  levels  of  damage. 


ii 


Figure  6.3  Cratering  and  penetration  events  in  Alnoo  targets  of  three  different 
thickness  (a=12.5  mm,  b=7.61  mm  and  c=l  .59  mm)  using  3.175  mm  diameter  soda-lime 
glass  projectiles  at  6  km/s  (Horz  et  al.  1995). 
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Table  6.1  Dimensional  characteristics  of  hypervelocity  impact  cases. 


Experimental 

PROJECTILE 

TARGET 

Case 

Diameter  (mm) 

Velocity  (km/s) 

Thickness  (mm) 

A 

3.175 

6 

12.5 

B 

3.175 

6 

7.61 

C 

3.175 

6 

1.59 

6.2.1  Computer  Simulations  Description 

The  simulations  in  Autodyn™  of  the  hypervelocity  impact  events  previously 
described  were  setup  using  Lagrangian  type  grids  for  both  the  target  and  projectile.  In  all 
three  cases  the  “Transmit”  boundary  condition  was  used  to  simulate  a  semi-infinite  plate 
widthwise.  The  “Transmit”  boundary  condition  allows  a  shock  wave  to  move  though  the 
boundary  of  the  grid  without  producing  any  reflection.  In  addition,  advantage  was  taken 
of  the  fact  that  the  impacts  were  normal  to  the  target  plate  allowing  axial  symmetry  to  be 
used  in  the  space  discretization  of  the  projectile  and  target.  Figures  6.4  to  6.6  depict  the 
grid  configuration  useu  for  the  three  simulations  corresponding  to  the  experiments 
previously  mentioned  (see  Figure  6.3  and  Table  6.1).  In  all  three  simulations,  the 
incremental-geometric-strain  “eroding”  algorithm  was  implemented  for  both  the  target 
and  projectile  with  a  2.5  limit  strain.  This  algorithm  “erodes”  (discards)  the  cells  that 
reach  the  limit  strain  thus  eliminating  the  cell-entangling  problem  that  exists  in 
Lagrangian  grids. 

Once  void  volume  fraction  reaches  the  critical  value  ^  (integration  of  Equation 
4.70)  at  any  cell  within  the  grid,  that  cell  is  “eroded”  away  to  simulate  the  local  fracture 
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of  the  material.  This  has  the  disadvantage  that  the  internal  energy  still  stored  in  the 
eroded  cell  is  eliminated.  Autodyn™  does  not  provide  any  other  way  to  simulate  fracture 
in  Lagrangian  grids,  and  the  implications  of  eroding  a  cell  with  high  internal  energy 
wannot  be  determined.  The  material  parameters  for  Alnoo  were  obtained  from  numerous 
published  works  of  other  authors.  Table  6.2  shows  the  material  parameters  for  Alnoo 
together  with  a  brief  description  of  where  they  appear  in  the  microdamage  constitutive 
model. 


Figure  6.4  Grid  configuration  for  simulation  (a)  corresponding  to  a  target  plate 
thickness  of  12.5mm. 


Figure  6.5  Grid  configuration  for  simulation  (b)  corresponding  to  a  target  plate 
thickness  of  7.61mm. 


Figure  6.6  Grid  configuration  for  simulation  (c)  corresponding  to  a  target  plate 
thickness  of  1.59mm. 
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Table  6.2  Material  parameters  for  Alnoo- 


Parameter 

Value 

Note 

Po 

2.77  g/cm3 

Initial  density 

Co 

5.328  km/s 

Bulk  sound  speed 

S 

1.338 

Slope  of  the  shock  to  particle  speeds 

To 

2.0 

Gruneisen  constant  j 

Cv 

875  J/Kg-°K 

Specific  heat 

To 

300  °K 

Initial  temperature  i 

|  T  mO 

923  °K 

Initial  melt  temperature 

k 

1.38x10'"  J/°K 

Boltzman  constant 

«  W 

0.9 

Fraction  of  plastic  work  dissipated  as  heat 

K0 

78.6  GPa 

Initial  bulk  modulus 

Po 

27.6  GPa 

Initial  shear  modulus 

A 

-7.245 

Constant  in  the  definition  of  the  Gruneisen  coefficient 

C. 

150  MPa 

Yield  stress  in  Johnson-Cook  (JC)  model 

i 

170  MPa 

Strain  hardening  constant  in  the  JC  model 

1  C3 

0.015 

Strain  rate  hardening  constant  in  the  JC  model 

n 

0.34 

Strain  hardening  exponent  in  the  JC  model 

S 

t  1.03 

Thermal  softening  constant 

i  a> 

1  1-5 

Constant  in  the  definition  of  the  melt  temperature 

a: 

0.004 

Elastic  coefficient  temperature  constant 

no 

j  10"  Pa-s 

Initial  material  viscosity 

Pn-.  |2xlO'JPa-s 

Melt  temperature  viscosity  i 

mi  j  2 

Exponent  in  the  plastic  rate  of  deformation  j 

nr3 

Initial  void  volume  fraction  j 

lr  >  0.3 

Critical  void  volume  fraction 

Maximum  void  volume  fraction  from  nucleation 

ao  j  10"*  cm 

Average  initial  radius  of  microvoids 

!  h 

7.54x10“  s" 

Constant  in  the  void  nucleation  model 

^no 

-300  Mpa 

Mean  stress  threshold  for  nucleation 

<p 

20 

Constant  in  the  void  growth  model 

m2 

10.3xl0'23 

Constant  in  the  void  nucleation  model 
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For  the  case  of  the  Soda-Lime  glass  projectile,  the  microdamage  constitutive 
model  presented  above  cannot  be  applied  since  the  model  is  applicable  only  for  ductile 
polycrystalline  materials.  For  this  reason,  the  Johnson-Holmquist  brittle  material  model 
already  implemented  in  Autodyn  was  used  with  the  materi"'  parameters  available  in  the 
materials  library  for  glass,  and  in  the  paper  by  Taylor  et  ul.  (1999). 

6.2.2  Simulation  Results 

The  graphics  capability  of  Autodyn™  allows  the  generation  of  contour  plots  of 
any  field  variable  that  forms  part  of  the  constitutive  microdamage  model.  For  the 
purpose  of  presenting  the  results  of  the  simulations  only  pressure,  temperature  and  void 
volume  fraction  contour  plots  will  be  shown. 

Figures  6.8  to  6.10  give  contour  plots  of  the  pressure,  temperature  and  void 
volume  fraction  for  the  12.5mm  thick  target  plate  (Case  A,  see  Table  6.1)  at  Ops,  0.5ps, 
l  .Ops,  1.5ps.  2ps,  2.5ps,  5ps  and  20ps  after  impact.  Figures  6.8(a),  6.9(a)  and  6.10(a) 
illustrate  the  initial  conditions  just  before  impact.  At  O.jps  after  impact  (see  Figures 
6.8(b).  6.9(b)  and  6.10(b)),  the  maximum  pressure  at  the  bottom  of  the  target  is  26GPa 
and  the  corresponding  temperature  is  680°K.  It  should  be  noted  that  as  the  shock- 
compression  wave  moves  through  the  target,  the  void  volume  fraction  is  being  reduced  to 
zero,  which  is  indicated  by  the  dark  blue  color  in  Figure  (6.10b).  At  the  moment  of 
impact,  the  pressure  in  the  target  reached  a  maximum  of  58GPa  and  the  corresponding 
temperature  is  900°K,  but  after  0.5ps  they  have  decayed  to  the  values  shown  in  the 
figures.  In  these  figures,  the  formation  of  a  target  rim  is  already  evident.  Figures  6.8(c- 
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e),  6.9(c-e)  and  6.10(c-e),  the  shock  wave  continues  to  move  across  the  target  with  the 
corresponding  decay  of  the  pressure  and  temperature  and  the  reduction  of  the  initial  void 
volume  fraction.  In  the  pressure  plots  (see  Figure  6.8),  it  should  be  noted  that  a  tensile 
wave,  reflected  from  the  rear  of  the  projectile  trails  the  shock-compression  wave.  Cnee 
the  shock  compression  wave  reaches  the  back  of  the  target  plate  it  is  reflected  as  a  tensile 
wave  that  is  superimposed  on  the  trailing  tensile  wave.  This  superposition  of  tensile 
waves  enhances  the  nucleation  and  growth  of  microvoids,  which  eventually  leads  to 
fracture.  This  effect  is  shown  in  Figures  6.8(e-f)  and  6.10(f),  where  the  void  volume 
fraction  has  reached  the  critical  value  and  spallation  (fracture)  has  begun.  Figures  6.8(h), 
6.9(h)  and  6.10(h)  illustrate  the  state  of  the  target  plate  at  the  end  of  the  simulation  (20ps 
after  impact).  At  this  time  in  the  simulation,  the  geometrical  changes  of  the  target  plate 
are  negligible  and  the  change  in  void  volume  fraction  is  zero.  The  predicted  Final  crater 
penetration  is  7.2mm  compared  to  7.3mm  in  the  experiment.  The  predicted  crater 
diameter  is  12mm  while  the  experimental  is  12.5mm.  This  gives  a  depth  of  penetration 
to  crater  diameter  (P/Dp)  of  0.6  for  the  simulation  and  0.58  for  the  laboratory  experment. 
The  predicted  diameter  of  the  spallation  area  is  12.5mm  and  is  19mm  for  the 
corresponding  laboratory  experiment,  indicating  a  thirty  three  percent  error  in  the 
prediction  of  the  spallation  diameter. 

Figures  6.11  to  6.13  give  contour  plots  of  the  pressure,  temperature  and  void 
volume  fraction  for  the  7.61mm  thick  target  plate  (Case  B,  see  Table  6.1)  at  Ops.  0.5ps. 

1  .Ops,  1 ,5ps,  2ps,  2.5ps,  5ps  and  20ps  after  impact.  Figures  6. 1 1  (a),  6. 1 2(a)  and  6. 1 3(a) 
illustrate  the  initial  conditions  just  before  impact.  The  evolution  of  the  simulations  for 
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Cases  A  and  B  (see  Table  6.1)  are  identical  up  to  approximately  l.Ops  after  impact, 
corresponding  to  contour  plot  (c)  of  Figures  6.8  to  6.10  and  6.11  to  6.13.  Figure  6.1 1(d) 
shows  the  initiation  of  spallation  after  the  reflection  of  the  compression  shock  wave  and 
its  superposition  with  the  trailing  tensile  wave.  At  2ps  after  impact  (Contour  plot  (e)  of 
Figures  6.1 1  to  6.13),  the  projectile  has  “eroded”  enough  of  the  target  plate  so  that  the 
crater  reaches  the  spallation  area.  In  contour  plots  (f)  to  (g),  the  spallation  diameter 
reaches  its  maximum  value  and  the  target  back  wall  continues  to  bulge,  driven  by  its 
momentum,  while  the  penetration  hole  gets  wider.  Figures  6.11(h),  6.12(h)  and  6.13(h) 
illustrate  the  state  of  the  target  plate  at  the  end  of  the  simulation  (20ps  after  impact).  The 
predicted  diameter  of  the  penetration  hole  is  1 1mm  while  for  the  laboratory  experiment  is 
1 1  4mm.  The  most  significant  difference  between  the  simulation  and  the  laboratory 
experiment  is  that  in  the  simulation  the  spall  does  not  detach  at  the  same  point  as  the 
laboratory  experiment.  The  fact  that  in  the  simulation  the  spall  section  that  is  left 
attached  to  the  target  plate  has  reached  an  equivalent  plastic  strain  of  1.2  (120%). 
indicates  the  possibility  that  another  fracture  mechanism,  one  that  is  not  taken  into 
account  in  the  simulation  and  is  operating  in  the  laboratory  experiment. 

Figures  6.14  to  6.16  give  contour  plots  of  the  pressure,  temperature  and  void 
volume  fraction  for  the  1.59mm  thick  target  plate  (Case  C,  see  Table  6.1)  at  Ops,  0.5ps, 
l.Ops,  1 ,5ps,  2ps,  2.5ps,  5ps  and  20ps  after  impact.  Figures  6.14(a),  6.15(a)  and  6.16(a) 
illustrate  the  initial  conditions  just  before  impact.  In  this  last  case,  the  projectile 
penetrates  the  target  plate  at  approximately  l.Ops  as  shown  in  contour  plot  (c)  of  Figures 
6.14  to  6.16.  After  initial  penetration,  the  penetration-hole  continues  to  grow  and  the 
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rims  at  the  front  and  back  of  the  target  plate  continue  to  develop  as  shown  in  contour 
plots  (d)  to  (g).  Figures  6.14(h),  6.15(h)  and  6.16(h)  illustrate  the  state  of  the  target  plate 
at  the  end  of  the  simulation  (20ps  after  impact).  The  predicted  diameter  of  the 
penetr.  ion.hole  is  8.5mm  while  for  the  actual  laboratory  experiment  it  is  8.8mm.  It  is 
interesting  to  note  the  microdamage  developed  close  to  the  wall  of  the  penetration  hole 
which  is  apparent  in  contour  plots  (c)  and  (d).  It  is  the  opinion  of  the  author  that  this  does 
not  happen  in  the  actual  laboratory  experiment,  but  it  is  not  possible  to  determine  if  this  is 
true  since  there  are  no  radiographs  at  different  times  of  the  target  plate  through  the 
penetration  process.  The  reason  for  the  appearance  of  microdamage  in  the  numerical 
simulation  mav  possibly  be  that  there  is  another  fracture  process  acting  that  is  not  being 
modeled  in  the  simulation.  For  example,  if  the  target  plate  fails  because  of  shear 
banding,  the  tensile  mean  stresses  needed  to  create  the  fracture  showTi  in  contour  plots  (c) 
and  (d)  may  not  be  reached.  This  phenomenon  is  demonstrated  by  Figure  6.7,  which 
shows  a  polished  and  etched  cross  section  through  a  rod  and  a  target  plate  for  a  case  of 
partial  penetration.  In  this  picture,  two  fracture  mechanisms  are  evident  in  the  target 
plate.  In  the  upper  portion  of  the  penetration  hole,  as  indicated  by  the  arrows,  shear 
banding  that  runs  through  the  remaining  section  of  the  plate  has  occurred,  while  in  the 
lower  section  of  the  hole,  shear  banding  is  not  evident  and  a  tensile  crack  has  developed 
instead.  Thus,  the  fact  that  formation  of  shear  bands  that  run  through  the  target  plate 
(plugging)  impedes  the  formation  of  cracks  or  microdamage  in  the  plate,  indicates  that 
the  implementation  of  a  shear  band  fracture  mechanism  is  needed  to  accurately  simulate 
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the  penetration  process  and  eliminate  the  unrealistic  development  of  microdamage  close 
to  the  wall  of  the  penetration  hole. 


Figure  6.7  Microfailure  features  observed  on  a  section  taken  through  a  steel  plate 
and  embedded  steel  rod  after  impact  (Shockey  et  al.  1979). 
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Figure  6.8(2)  Pressure  contour  plots  for  12.5mm  thick  Alnoo  target  plate  at  2ps, 
5ps.  5ps  and  20ps  (e,  f,  g,  h  respectively)  after  impact. 
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Figure  6.9(1)  Temperature  contour  plots  for  12.5mm  thick  Alnoo  target  plate  at 
Ops.  O.Sps,  l.Ops,  1.5ps  (a,  b,  c,  d  respectively)  after  impact. 
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Figure  6.9(2)  Temperature  contour  plots  for  12.5mm  thick  Alnoo  target  plate  at 
2ps.  2.5us.  5ps  and  20ps  (e,  f,  g,  h  respectively)  after  impact. 
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Figure  6. 10(1)  Void- volume  fraction  contour  plots  for  12.5mm  thick  Alnoo  target 
plate  at  Ops.  0.5ps,  l.Ops,  1.5ps  (a.  b,  c,  d  respectively)  after  impact. 
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Figure  6.1 1(1)  Pressure  contour  plots  for  7.61mm  thick  Alnoo  target  plate  at  Ops 
MSl  l.Ops.  1.5|is  (a,  b,  c,  d  respectively)  after  impact. 
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Figure  6.1 1(2)  Pressure  contour  plots  for  7.61mm  thick  Aluoo  target  plate  at  2ps, 
2.5ps,  5ps  and  20ps  (e,  f,  g,  h  respectively)  after  impact. 
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Figure  6.12(1)  Temperature  contour  plots  for  7.61mm  thick  Alnoo  target  plate  at 
Ops,  0.5ps,  l.Ops,  1.5ps  (a,  b,  c,  d  respectively)  after  impact. 
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Figure  6.12(2)  Temperature  contour  plots  for  7.61mm  thick  Alnoo  target  plate  at 
2 ps,  2.5ps,  5ps  and  20(is  (e,  f,  g,  h  respectively)  after  impact. 
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Figure  6.13(1)  Void-volume  fraction  contour  plots  for  7.61mm  thick  Alnoo  target 
plate  at  Ops,  0.5ps,  l.Ops,  1.5ps  (a,  b,  c,  d  respectively)  after  impact. 
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(e)  (f) 


(g)  (h) 


Figure  6.13(2)  Void-volume  fraction  contour  plots  for  7.61mm  thick  Alnoo  target 
plate  at  2|is.  2.5|is,  5ps  and  20ps  (e,  ft  g,  h  respectively)  after  impact. 
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Figure  6.14(1 )  Pressure  contour  plots  for  1.59mm  thick  Alnoo  target  plate  at  Ops, 
0.5gs,  l.Ops,  1.5ps  (a,  b,  c,  d  respectively)  after  impact. 
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Figure  6.15(1)  Temperature  contour  plots  for  1.59mm  thick  Aluoo  target  plate  at 
Ops,  0.5  ps,  l.Ops,  1.5ps  (a,  b,  c,  d  respectively)  after  impact. 
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Figure  6.15(2)  Temperature  contour  plots  for  1.59mm  thick  Al,ioo  target  plate  at 
,5ps,  5ps  and  20ps  (e,  f,  g,  h  respectively)  after  impact. 
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Figure  6.16(1)  Void- volume  fraction  contour  plots  for  1.59mm  thick  Aluoo  target 
plate  at  Ops.  0.5ps,  l.Ops,  1.5ps  (a,  b,  c,  d  respectively)  after  impact. 
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Figure  6.16(2)  Void-volume  fraction  contour  plots  for  1.59mm  thick  Alnoo  target 
plate  at  2ps,  2.5ps,  5ps  and  20ps  (e,  f,  g,  h  respectively)  after  impact. 


CHAPTER  7 


CONCLUSIONS  AND  RECOMMENDATIONS 

7.1  Introduction 

The  literature  review  presented  in  Chapter  1  demonstrated  the  necessity  to 
improve  upon  the  constitutive  models  used  for  the  simulation  of  high  shock  loading 
problems,  especially  in  relation  to  hypervelocity  impact.  As  was  stated  earlier,  the 
development  of  advanced  constitutive  models  for  the  simulation  of  hypervelocity  impacts 
is  of  primary'  importance  to  the  Aerospace  and  Defense  industry.  For  this  reason,  the 
goal  of  the  research  work  presented  in  this  dissertation  was  to  develop  a  constitutive- 
microdamage  model  that  can  describe  thermo-mechanical  deformation  associated  with 
high  shock  compression,  release,  dilatation  and  shear,  and  the  progressive  microdamage 
that  leads  to  fracture  of  ductile  polycrystalline  materials. 

7.2  Summary 

At  the  beginning  of  this  dissertation,  a  review  of  appropriate  literature  was 
conducted  and  a  brief  description  of  the  most  relevant  papers  was  presented  in  Chapter  1 . 

The  first  stage  in  the  development  of  the  constitutive  model  was  the  formulation 
of  the  Rankine-Hugoniot  relations  presented  in  Chapter  2,  which  were  derived  from  the 
“jump  conditions”  across  a  shock  wave  and  coupled  with  the  Us-Up  (shock  speed  to 
particle  speed)  empirical  linear  relation.  The  Rankine-Hugoniot  relations  were  used  to 
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give  the  final  form  of  the  Mie-Gruneisen  Equation  of  State,  which  was  developed  with 
detail  in  Chapter  3. 

The  focal  point  of  the  work  presented  in  this  dissertation  is  the  development  of  the 
microdamage  constitutive  model  that  describes  the  thermo-elastic  and  thermo¬ 
viscoplastic  behavior  of  polycrystalline  materials,  as  well  as  microdamage  evolution  in 
the  form  of  microvoid  nucleation  growth  and  coalescence.  The  development  of  the 
microdamage  constitutive  models  is  presented  in  Chapter  4. 

Due  to  the  limitations  imposed  by  the  numerical  framework  upon  which  the 
Autodyn  code  is  built,  a  series  of  simplifications  had  to  be  imposed  on  the 
implementation  of  the  microdamage  constitutive  model.  These  simplifications  affected 
the  plastic  rate  of  deformation  equation,  the  microvoid  growth  rate  equation  and  the 
temperature  rate  equation  in  the  manner  described  in  Chapter  5. 

To  test  the  performance  of  the  microdamage  constitutive  model,  a  series  of  three 
hvpervelocity  impact  laboratory  experiments  were  selected  for  simulation.  These  three 
experiments  covered  a  wide  range  of  projectile  diameter-target  thickness  rations.  In  the 
three  impact  simulations  Altoo  target  plates  were  impacted  with  soda-lime  glass  spherical 
projectiles  at  6.0  km/s.  A  detailed  description  of  the  laboratory'  experiments,  simulation 
parameters  and  simulation  results  were  given  in  Chapter  6. 
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7.3  Conclusions 

The  following  conclusions  can  be  reached  from  the  research  work  discussed  in  the 
previous  chapters: 

a)  The  microdamage  constitutive  model  as  implemented  in  the  Autodyn™  code 
successfully  simulated  the  depth  of  penetration  and  crater  diameter  for  the 
1 2.5mm  thick  target  plate  (Case  A)  with  an  error  below  five  percent. 

b)  The  diameter  of  the  spall  for  the  12.5mm  thick  target  plate  (Case  A)  was 
under  predicted  by  thirty  three  percent.  It  is  the  view  of  the  author  that  this 
error  could  be  reduced  if  some  of  the  material  parameters  presented  in  Table 
6.2  were  modified.  However,  the  decision  was  made  to  keep  the  material 
parameters  within  the  ranges  listed  in  the  literature. 

c)  The  constitutive  model  successfully  predicted  the  diameter  of  the  penetration 
hole  in  the  7.61mm  thick  plate  (Case  B)  to  within  a  four  percent  error. 

d)  There  is  a  discrepancy  in  the  extent  of  spall  detachment  for  the  7.61mm  thick 
plate  (Case  B)  between  numerical  simulation  and  the  laboratory 
experiment,  but  in  general,  the  main  deformation  and  fracture  features  are 
reasonably  simulated. 

e)  The  constitutive  model  successfully  predicted  the  diameter  of  the  penetration 
hole  in  the  1.59mm  thick  plate  (Case  C)  to  within  a  four  percent  error. 

f)  The  existence  of  microdamage  developed  close  to  the  wall  of  the  penetration 
hole  in  the  1.59mm  thick  plate  (Case  C)  cannot  be  corroborated,  since  there 
are  no  radiographs  available  of  the  state  of  the  target  plate  at  approximately 
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1 .Ops  after  impact  for  the  actual  laboratory  experiment.  Even  though  it  is  not 
possible  to  corroborate  the  existence  of  microvoid  damage,  the  fact  that  the 
numerical  simulation  shows  rough  surfaces  on  the  wall  of  the  penetration 
hole,  while  the  experiment  does  not,  seems  to  indicate  that  the  micro  damage 
process  might  not  be  present  in  the  actual  laboratory  experiment, 
g)  The  discrepancies  between  the  numerical  simulations  and  experimental  results 
described  in  conclusions  (d)  and  (f)  seem  to  indicate  that  there  is  the 
possibility  that  another  failure  mechanism,  such  as  shear  band  formation;  that 
is  not  taken  into  account  in  the  model  simulation,  is  present  in  the  laboratory 
experiment. 

7.4  Recommendations 

Because  of  the  discrepancies  mentioned  previously,  it  is  recommended  that  the 
microdamage  evolution  model  be  complemented  by  a  shear  or  shear  band  failure  model. 
This  would  improve  the  constitutive  model  presented  here  and  would  r.uend  its 
successful  applicability  to  a  wider  range  of  projectile  diameter  to  target  plate  thickness 
ratios. 

It  is  also  recommended  that  the  microdamage  constitutive  model  be  implemented 
for  a  Eulerian  type  grid.  This  would  eliminate  the  necessity  to  implement  an  “eroding” 
algorithm  to  eliminate  highly  deformed  cells,  and  to  simulate  fracture  in  the  form 
described  in  Chapter  6.  Furthermore,  the  possibility  of  doing  simulations  in  a  three- 
dimensional  configuration  should  be  considered.  This  should  be  done  because  in  the 
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actual  laboratory  experiment,  the  axial  symmetry  is  only  present  until  the  spall  detaches 
and  petal  shape  fracture  begins  to  develop.  This  could  be  the  reason  for  the  discrepancies 

in  the  detached  spall  fracture  mentioned  in  conclusion  (d). 

The  microvoid  nucleation  model  is  considered  by  the  author  to  be  one  of  the  least 
developed  components  of  the  microdamage  constitutive  model  presented  here.  In  the 
literature  review  conducted  as  part  of  the  research  program,  there  are  no  reports  of  studies 
of  the  nucleation  process  under  mean  stresses  that  are  at  least  one  order  of  magnitude 
greater  that  the  nucleation  threshold,  and  at  strain  rates  of  the  order  of  108  s  ',  as  is  the 
case  in  the  hvpervelocity  impact  experiments  discussed  in  the  simulations  part  of  this 
dissertation.  Even  though  the  microvoid  nucleation  model  was  improved  here  by 
incorporating  a  strain  and  strain  rate  dependent  nucleation  threshold,  it  is  the  view  of  the 
author  that  further  research  has  to  be  conducted  on  this  subject. 

As  already  mentioned,  the  energy  rate  equation  (4.86)  was  derived  trough  the  use 
of  the  energy  rate  balance  equation  that  applies  to  any  material  under  any  loading 
condition,  but  it  presented  difficulties  in  its  numerical  implementation  as  described  in 
Section  5.1.3,  and  the  model  already  implemented  in  Autodyn™  had  to  be  used  instead. 
Furthermore,  it  was  demonstrated  in  Section  6.4  that  there  is  a  significant  difference  in 
the  temperatures  calculated  with  the  energy  rate  balance  equation  (4.86)  and  the 
calculation  with  classical  thermodynamics,  upon  which  the  Autodyn™  model  is  based  on 
in  part.  For  this  reason,  further  research  has  to  be  conducted  to  determine  which 
calculation  predicts  the  actual  experimental  temperature  under  hypervelocity  impact  more 
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accurately.  This  is  a  challenging  problem  because  of  the  difficulties  in  measuring 
temperature  fluctuations,  in  laboratory  experiments,  that  last  only  a  few  nanoseconds. 

Lastly,  it  is  recommended  that  the  model  be  tested  for  different  projectile  and 
target  materials.  This  would  give  more  confidence  in  the  implementation  of  the 
microdamage  constitutive  model  presented  in  this  dissertation. 
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